ANALYSIS AND INVERSION OF SEISMIC REFRACTION TRAVELTIMES by DAVID FRANKLIN ALDRIDGE B.A., Amherst College, 1975 M.A., The University of California, Berkeley, 1980 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES Department of Geophysics and Astronomy We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April 1992 Β© David Franklin Aldridge, 1992 Signature(s) removed to protect privacy 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. (Signature) __________________________ Department of wc/ A%cQ4o,.1 7 / The University of British Columbia Vancouver, Canada Date 27 DE-6 (2/88) Signature(s) removed to protect privacy 11 ABSTRACT Techniques for forward modeling and inversion of head wave traveltimes within the framework of one and two dimensional earth models are well developed. The first portion of this thesis extends these methods to encompass three dimensional layered models. Each critically refracting horizon of the model is approximated by a plane interface with arbitrary strike and dip. An advantage of this simple representation is that rapid computation of head wave traveltimes for arbitrary source-receiver geometries can be achieved with a minimum of ray tracing. Inversion methods are then developed for estimating the parameters defining single-layer and multilayer earth models. For the single-layer model, an algebraic solution to the inverse problem exists if refraction traveltimes are observed along two independent line profiles. For multilayer models and/or nonproffle recording geometries, the inversion is formulated as a constrained parameter optimization problem and solved via linear program ming. Inclusion of constraints, in the form of inequality relations satisfied by the model parameters, often governs the ability of the algorithm to converge to a realistic solution. The procedure is tested with traveltimes recorded on broadside profiles in a deep crustal seismic experiment. The second part of this thesis provides specific improvements to various two dimen sional refraction traveltime inversion techniques. The generalized reciprocal method (GRM) is reformulated on the basis of an earth model characterized by vertical, rather than nor mal, layer thicknesses. This allows point values of interface depth to be inferred from the observed traveltimes. A novel interpretation method (critical offset refraction profiling) is described that yields point values of interface depth, interface dip, and refractor velocity. A smooth depth profile of the refracting horizon is then constructed using techniques of linear inverse theory. Finally, an automated version of the classical wavefront method for inter preting refraction traveltimes is developed. Recorded arrival times are downward continued through a near surface heterogeneous velocity structure with a finite-difference propagation algorithm. The locus of a refracting horizon is then obtained by applying a simple imaging condition involving the reciprocal time (the source-to-source traveltime). The method is UI tested, apparently successfully, on a shallow refraction dataset recorded at an archeological site. The final portion of this thesis develops an iterative tomographic inversion procedure for reconstructing a two dimensional P-wave velocity field from measured first arrival times. Two key features of this technique are (i) use of a finite-difference algorithm for rapid and ac curate forward modeling of traveltimes, and (ii) incorporation of constraint information into the inversion to restrict the nonuniqueness inherent in large scale, nonlinear tomographic inverse problems. Analysis of a simulated vertical seismic profile (VSP) plus crosswell exper iment indicates that the inversion algorithm can accurately reconstruct a smoothly varying interwell velocity field. Inclusion of constraints, in the form of horizontal and vertical first- difference regularization, allows the solution of a traveltime tomography problem that is otherwise severely underdetermined. iv TABLE OF CONTENTS Abstract . ii Table of contents iv List of tables . . . . viii List of figures 1 Introduction 2 Head wave traveltimes in a three dimensional single-layer earth 2.1 Introduction 2.2 Earth model and recording geometry 2.3 Traveltime derivation 2.4 Traveltime inversion 2.4.1 Special cases 2.4.2 Arbitrary line azimuths 2.5 Traveltime isochrons 2.6 Inversion of traveltime isochrons 2.7 Conclusion 3 Head wave traveltimes in a three dimensional multilayered earth 3.1 Introduction 3.2 Earth model 3.3 Raypath geometry 3.4 Traveltime derivation 3.4.1 Downgoing traveltime 3.4.2 Upgoing traveltime 3.4.3 Critically refracted traveltime 3.4.4 Total traveltime 3.4.5 Variants of the basic formula 10 10 10 β’ . . 13 17 18 20 β’ . . 21 26 27 29 29 30 β’ β’ . 32 36 β’ . . 36 39 β’ . 40 β’ . 40 β’ . . 42 3.5 Generalizations of the traveltime formula . 3.5.1 Source and receiver on separate interfaces 3.5.2 Arbitrary source and receiver locations 3.5.3 Arbitrary reference points for layer thickness 3.5.4 Reflection traveltime 3.6 Rapid traveltime computation 3.6.1 General description 3.6.2 Raytracing technique 3.7 Forward modeling examples 3.7.1 Profile geometry 3.7.2 VSP geometry 3.8 Conclusion V 4444464852545557 61 61 64 66 4 Inversion of head wave traveltimes for three dimensional planar structure 4.1 Introduction 4.2 Inversion mathematics 4.2.1 General theory 4.2.2 Calculation of sensitivities . 4.3 Synthetic examples 4.3.1 Single layer 4.3.2 Multiple layers 4.4 Field data example 4.4.1 Peace River Arch broadside data 4.4.2 Static corrections 4.4.3 Inversion results 4.5 Conclusion 68 6869697276 76 80818183 87 94 5 Point depth and dip estimates from refraction traveltimes 5.1 Introduction 96 96 5.2 Specialization to a 2D model 5.3 Derivation of GRM parameters 5.3.1 Velocity analysis function 5.3.2 Apparent horizontal velocity 5.3.3 Generalized time-depth 5.3.4 Depth conversion factor 5.3.5 Time-depth near a shotpoint 5.4 Critical offset refraction profiling 5.4.1 Earth model 5.4.2 Inversion method 5.4.3 Critical offset determination 5.5 Conclusion vi 97 101 101 103 108 110 112 113 113 115 118 120 7 Refractor imaging using an automated wavefront reconstruction method 7.1 Introduction 7.2 Finite-difference traveltimes 7.2.1 Wavefront construction 122 122 124 125 127 128 130 133 134 137 144 6 Construction of a smooth refractor depth profile 6.1 Introduction 6.2 Smoothest model construction 6.2.1 Modified data equation 6.2.2 Objective function 6.2.3 Extremizing the objective function . . 6.2.4 Constructing the model 6.3 Interpolation via the smoothest model 6.3.1 General theory 6.3.2 Numerical example 6.4 Conclusion 146 146 147 147 vu 7.2.2 Wavefront reconstruction 7.3 Refractor imaging 7.4 Refractor velocity estimation 7.5 Field data example 7.6 Conclusion 152 152 160 162 167 8 Two dimensional tomographic inversion with finite-difference traveltimes 8.1 Introduction 8.2 General theory and model representation 8.3 Forward modeling 8.4 Ray generation 8.5 Inversion mathematics . 8.6 Model constraint equations 8.7 Synthetic examples 8.8 Conclusion 169 9 Summary References Appendices A Conic sections in polar coordinates B Traveltime analysis for a simple 3D model C Traveltime analysis for a simple 2D model D Derivation of formula (6.13) E Vanishing of two derivatives F Spatially correlated traveltime errors . G Ray initiation and termination logic . H Optimum starting slowness 193 196 202 202 205 208 212 213 216 218 221 169 171 172 175 178 182 184 189 VIII LIST OF TABLES 8.1 Example of ray tracing logic 177 Cl Ray initiation logic for a receiver located within a cell 220 ix 2.1 Single layer earth model and coordinate reference frame . 2.2 Plan view of surface recording geometry 2.3 Head wave raypath 2.4 Normalized apparent refractor velocity 2.5 Head wave traveltime isochrons 3.1 Orientation of i interface of multilayered model - 3.2 Schematic diagram of a raypath critically refracted at kth interface 3.3 Snellβs law of refraction at th interface 3.4 Snellβs law of reflection at kthl interface 3.5 Schematic representation of a critically reflected/refracted raypath 3.6 Direct and head wave traveltimes recorded by four reversed profiles 3.7 Head wave traveltimes recorded in a VSP configuration 4.1 Plan views of two areal recording arrays 4.2 Broadside recording geometry for Peace River Arch (PRA) experiment 4.3 Receiver static functions for PRA lines A and B 4.4 Comparison between PRA observed and predicted traveltimes 4.5 Plan view of vertical depths to Moho in PRA region 4.6 Vertical depths to Moho beneath PRA lines A and B 5.1 Critically refracted raypaths for ORM analysis 102 5.2 Comparison of ORM apparent velocities with true velocity 106 5.3 Two dimensional model used for critical offset refraction profiling 114 5.4 Expanded view of triangle ABC 117 5.5 Determination of critical offset distance 119 6.1 Refractor elevation model constructed from accurate depth samples 138 6.2 Refractor elevation model with optimized parameters d1 and d2 . . . . . 139 LIST OF FIGURES 11 12 13 16 24 31 33 35 53 56 62 65 77 82 88 90 91 94 6.3 Refractor elevation model constructed from accurate depth and dip samples 141 6.4 Refractor elevation model constructed from inaccurate depth and dip samples x 142 7.1 Subsurface first arrival wavefronts calculated with finite-differences 7.2 Surface arrival time curves 7.3 Finite-difference traveltime extrapolators 7.4 First arrival wavefronts over a shallow syncine 7.5 Reconstructed first arrival wavefronts over a shallow syndine 7.6 Imaging a syndinal refractor 7.7 Imaging an antidinal refractor 7.8 Dependence of syncine on imaging time and downward continuation 7.9 Reconstructed wavefronts and syncline image; noisy traveltimes 7.10 Refractor velocity estimates for syncline and anticine models 7.11 Shallow refraction traveltime data from Phalasarna, Crete . 7.12 Shallow refractor image from Phalasarna data 7.13 Comparision of observed and predicted traveltimes 7.14 Forward and reverse first arrival wavefronts for final earth model 148 149 151 153 154 156 157 velocity . . 158 160 161 163 164 165 166 8.1 Example raypath through cell ij of a 2D velocity model 8.2 A 2D velocity model with wavefronts and raypaths from a surface source 8.3 Five-cell constraint operator 8.4 243 double-well VSP and crosswell raypaths 8.5 Velocity tomograms obtained from VSP and crosswell arrivals; accurate data 8.6 Velocity tomogram obtained from VSP and crosswell arrivals; noisy data 8.7 81 crosswell raypaths 8.8 Velocity tomograms obtained from crosswell arrivals only Al A noncentered, rotated effipse Cl A simple 2D earth model with three layers Gi Grid cells surrounding a source or receiver point 176 179 β’ 183 β’ 186 β’ 187 β’ . 188 β’ . 190 β’ . 191 β’ . 204 β’ . 209 219 1CHAPTER 1 INTRODUCTION The traveltimes of head waves propagating in layered earth models have been studied extensively since the inception of applied seismology in the 1920βs. The earliest English language publications that examine this topic appear to be the classic works by Barton (1929) and Heiland (1929). These authors derive head wave traveltinie formulae appropriate for a simple model consisting of a layer overlying a haifspace. The two media have uniform - P-wave velocities and are separated by a plane interface. If the critically refracting interface is horizontal, then head wave traveltime is a linear function of source-receiver offset distance X: T(X) = mX + b. (1.1) The slope m and intercept b are independent of the recording geometry and are easily determined from the known earth model parameters. Numerous investigators, beginning with Barton (1929), have extended this eicpression to multilayered one dimensional models. If the subsurface interface is clipping, then the same general mathematical formula ap plies. However, the slope and intercept are no longer invariant. Rather, the intercept depends on the horizontal coordinates of the source (xs, ys) and the slope depends on the azimuth angle Iβ of the receiver relative to the source. Barton (1929) and Heiland (1929) examine the case where an inline array of receivers is oriented normal to the strike direction of the refracting horizon. Hence, the azimuth angle EI is restricted to two values (say, = 0 and βIβ = ir) that collectively refer to the updip and downdip directions. Surface-to-surface head wave traveltime is T(X,,x) = m(lβ)X + b(x). (1.2) In this two dimensional situation, the plane of the head wave raypath is vertical. Barton (1929) and Heiland (1929) derive closed form mathematical expressions for the slope and 2intercept in terms of the parameters defining the single-layer earth model. Thus, equation (1.2) is easily evaluated for this case. Head wave traveltime analysis is extended to multilayered two dimensional models with plane, dipping interfaces by Ewing et al. (1939), Dooley (1952), Adachi (1954), Mota (1954), Ocola (1972), Johnson (1976), and Diebold and Stoffa (1981). The linear traveltime expression (1.2) also applies in this situation. Once again, the recording profile must be oriented perpendicular to the common strike direction of the subsurface refracting horizons. If this condition is not maintained, then the raypath is not confined to a single vertical plane, and a full three dimensional treatment of the problem is necessary. For the multilayered case, the slope and intercept in (1.2) cannot be simply expressed in terms of the specified earth model parameters. Rather, each depends implicitly on the model parameters through a set of raypath orientation angles. It is possible to determine these angles via repeated application of Snellβs law of refraction, starting at the critically refracting horizon and working upward to the surface. Hence, equation (1.2) can be evaluated numerically. Methods for inverting observed head wave arrival times to obtain the earth model pa rameters are also discussed by the above authors. The techniques assume that traveltimes are recorded by inline forward and reverse spreads, or by a split spread. If the P-wave ve locity of the uppermost layer is known, then the measured slope of a particular traveltime branch can be used to infer the angle of incidence of the head wave at the surface. Repeated application of Snellβs law of refraction from the surface downward determines the raypath orientation angle within the layer immediately above the critically refracting interface. The dip and velocity of the critical refractor are then obtained from raypath orientation angles calculated for forward and reverse propagating head waves. The method assumes that all velocities and dip angles overlying the target horizon haye been previously determined. Fi nally, interface depths below the source locations are calculated from the measured intercept times. An alternative two dimensional inversion method exploits the variation of intercept time b(x) with source position. Cunningham (1974) demonstrates that head wave arrival times observed on two inline spreads with the same source-receiver azimuth βP can be analyzed to 3yield refractor dip, depth, and velocity. Although the slopes of the two traveltime curves are the same, the intercept times associated with the different source locations provide sufficient independent information for a solution. The method has practical utility in marine seismic exploration, where refracted arrivals are recorded on single-ended (i.e., unreversed) spreads. Cunningham (1974) examines simple models consisting of one and two layers overlying a halfspace. However, it is possible to demonstrate the validity of the method for general multilayered models. The fact that both forward and reverse arrivals are not required for a successful two dimensional inversion is not widely appreciated, even by specialists in seismic -refraction exploration (e.g., Lankston, 1990; Palmer, 1991). The purpose of Chapters 2, 3, and 4 of this thesis is to extend forward modeling and inversion of head wave arrival times to three dimensional layered models. This class of earth models is characterized by uniform velocity layers bounded by plane interfaces pos sessing arbitrary dip and strike. The linear traveltime expression remains valid in this three dimensional situation: T(X,,xs,ys) = m(P)X + b(xs,ys,β). (1.3) The source-receiver azimuth angle is not restricted to two values, but may range throughout the interval 0 <2ir. Note that the slope in equation (1.3) depends only on the recording azimuth prβ’ In general, the intercept depends on both the source location and the azimuth angle to the receiver. However, for the simple model consisting of a single layer overlying a halfspace, the dependence on βP vanishes. Chapter 2 presents a rigorous derivation of equation (1.3) for the single-layer earth model. Closed form mathematical expressions are obtained for the slope and intercept in terms of the specified model parameters. Also, equations for critical and crossover distances are derived. All of these expressions are generalizations of more- familiar formulae that apply to one and two dimensional models, and are quite useful for forward modeling purposes. Unlike the analogous situation for seismic reflection, there is very little published information on the traveltimes of critically refracted waves in this basic three dimensional model. Dix 4(1935), using geometric arguments, demonstrates that head wave isochrons (i.e., loci of fixed arrival time at the surface) are conic sections. A similar result is obtained by Dunkin and Levin (1971) although their analysis is limited to a particular areal recording configuration. Chapter 2 provides a derivation of the same result, using a completely different methodology from that in Dix (1935). No restrictions are placed on the recording geometry. The simple βlayer over a halfspaceβ earth model is defined by five parameters: two P-wave velocities, two interface orientation angles, and a depth to the interface below a reference point. If the overburden velocity vi is known, then only four model parameters need to be determined by inverting the observed head wave arrival times. Heiland (1940, p. 525-527) describes an inversion method that utilizes traveltimes measured by four line profiles ema nating from a single source location. Russell et al. (1982) require three proffles. However, Chapter 2 demonstrates that traveltimes recorded by only two proffles are sufficient to deter mine the attitude, velocity, and depth of the refractor. Two independent line profiles provide four measured data (two slopes and two intercepts) that allow an algebraic solution for the four unknown model parameters. Inversion methods that are not based on conventional line proffle recording geometries are also possible. For example, Chapter 2 indicates that all five model parameters can be obtained by analysis of the geometric properties of head wave isochrons. This result suggests interesting possibilities for determining overburden velocity v from the refraction data alone. In Chapter 3, the derivation of the linear traveltime formula (1.3) is extended to multi layered earth models. A substantial simplification in the mathematical proof is obtained by using a novel three dimensional form of Snellβs law of refraction. A previous proof given by Diebold (1987) is very cumbersome, and may actually be incorrect because it yields an ex pression that does not generalize to arbitrary source-receiver geometries properly. As in the two dimensional multilayered case, the slope and intercept in (1.3) depend on the orientation angles of the raypath within each layer. However, two angles are now necessary to describe the orientation of a raypath segment in three dimensional space. Chapter 3 describes a rapid computational procedure for obtaining these angles by applying Snellβs law repeatedly from 5the critically refracting horizon upward to the surface. Hence, the traveltime equation (1.3) is readily evaluated. Finally, inversion of head wave traveltimes for three dimensional planar structure is ad dressed in Chapter 4. Single layer and multiple layer algorithms are developed for data acquisition geometries where sources and receivers are located on the surface. The inver sion is posed as a constrained parameter optimization problem. An initial estimate of the earth model parameters is iteratively refined until an acceptable match is obtained between observed and predicted arrival times. Similar parametric inversion schemes have been re cently applied to three dimensional reflection traveltimes by Chiu et al. (1986), Lin (1989), and Phadke and Kanasewich (1990). Chiu and Stewart (1987) invert combined VSP and surface-to--surface reflection traveltimes. Kanasewich and Chiu (1985) invert combined re flection and refraction arrivals. They use the iterative ray-bending approach of Chander (1977a) to calculate the traveltime derivatives needed for the inversion. The computational procedure developed in Chapter 3 is not iterative. Hence, it provides a rapid alternative to ray-bending or ray-shooting methods for obtaining the necessary head wave traveltime derivatives. A novel feature of the inversion algorithms discussed in Chapter 4 is the introduction of constraint information via inequality relations that are imposed on the model parame ters. Often, a priori geological or geophysical data are available to guide and constrain an inversion. For example, many shallow seismic refraction projects are undertaken in conjuc tion with a drilling program. Interface depth and layer velocity data may be obtained from borehole logs. Constraints are particularly useful for the inversion of head wave traveltimes, because the problem can be very ill-conditioned and admit numerous solutions. Application of the constrained inversion algorithm to a crustal seismic refraction dataset from northern Alberta indicates that it can be a useful tool for analysis of broadside recorded traveltimes. Currently, there is a lack of techniques for the effective analysis of such data. Following the analysis of three dimensional earth models in Chapters 2, 3, and 4, the central portion of this thesis (Chapters 5, 6, and 7) develops specific improvements to various two dimensional refraction traveltime inversion methods. By restricting consideration to two 6dimensional models, examination of geological realities like undulating interfaces and/or variable velocity media becomes mathematically tractable. Chapter 5 demonstrates that point estimates of interface depth and dip can be inferred from the observed refraction traveltimes. Two interpretation procedures are used for the discussion. The Generalized Reciprocal Method (GRM) is a technique for delineating an undulating subsurface interface from refracted arrivals recorded on inline forward and reverse spreads. It was developed by Palmer (1980, 1981) as an extension of the conventional reciprocal method (Hawkins, 1961) of refraction traveltime interpretation. Although the GRM has been successfully applied to the problem of mapping undulating horizons, the mathematical formulation of the method is based on a two dimensional earth model with plane interfaces. Hence, the improved head wave traveltime formula derived in Chapter 3 can be applied to GRM analysis. A useful traveltime expression is obtained by specializing the general equation to a two dimensional situation. As demonstrated in Chapter 5, all of the common GRM analysis tools can be derived from this novel 2D traveltime formula in a straightforward manner. Moreover, these new expressions allow point depth estimates of the refracting interface to be calculated from the measured traveltimes. A refractor depth profile can then be obtained by interpolation. In contrast, the current GRM depth estimation procedure involves the laborious task of constructing an envelope to a set of circular arcs (Hatherly, 1980). A traveltime inversion technique that explicitly incorporates nonplane refracting hori zons into the model should yield more accurate results. Thus, Chapter 5 also proposes a head wave traveltime interpretation method tentatively named critical offset refraction profiling. This inversion technique accommodates undulations in both the surface and the refracting interface, as well as horizontal variations in the velocity of the refracting layer. Analysis reveals that a mathematically exact solution for the earth model parameters exists: point values of interface depth, interface dip, and refractor velocity can all be obtained from the observed traveltimes. Hence, as with the modified GRM developed in the same chapter, a continuous depth profile of the interface can be constructed by interpolation. The dip 7estimates obtained from critical offset refraction profiling provide additional constraints on the interpolant. A detailed description of a method for calculating a smooth refractor depth profile from a set of point estimates of depth and dip is presented in Chapter 6. Determination of the interface depth proffle is treated as an interpolation problem. Hence, the technique differs significantly from traditional methods of depth proffle computation via envelope construction (Thornburgh, 1930; Dix, 1941; Tarrant, 1956; Hales, 1958; Hawkins, 1961; Hatherly, 1980; Palmer, 1980, 1990, 1991). Linear inverse theory provides the mathematical framework for solving this problem. The smoothest (i.e., minimum curvature) interpolant is the natural model to adopt in this situation because refraction traveltime inversion methods assume, either explicitly or implicitly, that local interface curvature is negligible. Thus, the final model for the refracting horizon is consistent with prior assumptions used for inferring its depth and/or dip. Posing the problem as an interpolation issue also has several specific advantages regarding treatment of the data: (i) additional depth and dip values arising from a variety of geological, geophysical, or engineering techniques are readily incorporated into the model construction, (ii) variable weighting of all data is easily achieved, and (iii) an adjustable misfit to error contaminated data is possible. Recent advances in seismic refraction interpretation involve constructing a subsurface image directly from the observed first break waveforms (Clayton and McMechan, 1981; Hill, 1987). Chapter 7 describes an alternative imaging procedure that works only with the picked arrival times. The technique is an automated implementation of the classical Wavefront Method for interpreting refraction traveltimes (Thornburgh, 1930; Hagedoorn, 1959; Rock well, 1967). Modern finite-difference propagation algorithms are used to downward continue recorded refraction arrival times through a near-surface heterogeneous velocity structure. Two such subsurface traveltime fields need to be reconstructed from the arrivals recorded on a forward and reverse geophone spread. The locus of a shallow refracting horizon is then defined by a simple imaging condition involving the reciprocal time (the traveltime between source positions at either end of the spread). Refractor velocity is estimated in a subsequent step by calculating the directional derivative of the reconstructed subsurface wavefronts 8along the imaged interface. The principal limitation of the technique arises from imprecise knowledge of the overburden velocity distribution. This velocity information must be ob tained from uphole times, direct arrivals, shallow refractions, and borehole data. Analysis of synthetic data examples indicates that the technique can accurately image both syncinal and antidinal structures. The method is also tested, apparently successfully, on a shallow refraction dataset acquired at an archeological site in western Crete. Recently, tomographic techniques have been applied to the shallow seismic refraction problem (Hampson and Russell, 1984; De Amorim et al., 1987; Olsen, 1989; Scales et al., 1990). Traveltime tomography is a general and versatile method of velocity analysis. The final portion of this thesis (Chapter 8) presents an iterative tomographic inversion procedure for reconstructing a two dimensional velocity field from measured first arrival times. Two key features of this technique are (i) use of a finite-difference algorithm for rapid and accurate forward modeling of traveltimes, and (ii) incorporation of constraint information into the inversion in order to restrict the nonuniqueness inherent in large scale, nonlinear tomographic inverse problems. Finite-difference traveltime computation (Vidale, 1988) provides a useful alternative to conventional raytracing in tomography. All first arrival wave types (direct and refracted arrivals, head waves, diffractions) are handled with relative ease. Curved raypaths, needed for subsequent updating of the velocity model, are generated by following the steepest descent direction through a computed traveltime field from each receiver back to the source. Since arrival times are computed throughout a gridded two dimensional velocity field, very general recording geometries are easily accommodated. The main limitation of the method is that it is restricted to first arrivals. Hence, the algorithm developed in Chapter 8 cannot be applied to reflection tomography. Constraint information may arise from known values of the subsurface velocity (e.g., from outcrops or well logs) or the imposition of reasonable geological characteristics, like smoothness, on the constructed velocity model. In either case, the data equations are aug mented with additional linear constraint equations and a least squares solution of the entire 9system is obtained. In Chapter 8, the constrained inversion algorithm is applied to a sim ulated double-well VSP plus crosswell experiment. The results indicate that the algorithm can accurately reconstruct a smoothly varying interwell velocity field. Inclusion of constraint information, in the form of horizontal and vertical first-difference regularization, allows the solution of a traveltime tomography problem that is otherwise severely underdetermined (243 data and 10000 unknowns). For this example, no obvious improvement is obtained by the addition of borehole velocity constraints. 10 CHAPTER 2 HEAD WAVE TRAVELTIMES IN A THREE DIMENSIONAL SINGLE-LAYER EARTH 2.1 Introduction Many refraction seismologists believe that a minimum of two reversed profiles, preferen tiafly oriented at right angles to each other, are necessary to determine the attitude, velocity, and depth of a plane subsurface refractor (e.g., Heiland, 1940, p. 525-527). Russell et al. (1982) demonstrate that traveltimes recorded along three unreversed spreads can be lyzed to yield this same information. They state that similiar measurements made by only two such profiles cannot define the three dimensional attitude of the dipping horizon. This statement, however, is incorrect. Russell et al. (1982) do not fully utilize the information contained in the intercept times of the traveltime curves. One purpose of this chapter is to demonstrate that, in many cases, two refraction proffles are sufficient to define the three dimensional attitude, true velocity, and depth of a plane refractor. Generally speaking, the main condition required is that the two lines provide independent traveltime information about the subsurface. 2.2 Earth model and recording geometry The earth model consists of a single layer with P-wave speed v overlying a halfspace with P-wave speed V2. The plane interface separating the two media possesses, in general, a three dimensional dipping attitude. Let the zy plane of a right handed, rectangular coordinate system be coincident with the free surface of the layer; depth z is measured positive in the downward direction. If r = xi + yj + zk is the position vector of an arbitrary point on the subsurface interface, then the equation defining this dipping plane is rn = d. (2.1) 11 d is the perpendicular distance from the coordinate origin 0 to the plane and n is a unit vector pointing from 0 along the normal to the plane. Figure 2.1 illustrates that n is conveniently described in terms of polar and azimuthal angles q and 0: = (sinc cos0)i+(sinΓ§t sin0)j-1-(cosq)k. (2.2) Γ§ (0 < ir/2) is the dip angle of the interface and 0 (0 9 < 2ir) defines the upc]ip direction. If the +z and +y axes are taken to point toward geographic north and east, respectively, then the strike angle of the interface is 6? + ir/2 (modulo 2ir). y 0 x n z Fig. 2.1. Earth model and coordinate reference frame. 12 Figure 2.2a is a plan view of the surface recording geometry. The position vectors of the source S and receiver R are given by = Xsi + ysj, rfl Xfl1 + YRJβ’ (2.3a,b) The azimuth angle of the receiver relative to the source is a (0 a < 2ir). Azimuth is measured positive in the clockwise sense from geographic north. If the horizontal offset between source and receiver is denoted by X (X 0), then the receiver coordinates can be expressed in terms of the source coordinates as follows: XR = as+Xcosa, YR = Ys+XSina. (2.4a, b) Proffle recording geometry is defined by the condition that the azimuth angle a remains fixed for a set of receivers that record energy from a single source. a) b) x x R 0 0 x S y y Fig. 2.2. (a) Plan view of surface recording geometry. S and R denote a point source and a point receiver, respectively. (b) Spatial relationship of the two sources S1 and 52. 13 2.3 Traveltime derivation Consider the critically refracted raypath from surface source S to surface receiver R (Figure 2.3). This raypath is confined to a single plane referred to herein as the raypath plane. In general, the raypath plane is not a vertical plane (parallel to the z axis). Only in the particular case where the profile line is oriented directly updlip or downdip is the raypath plane vertical. The traveltime of a head wave propagating along the critical raypath from source to receiver can be worked out from simple geometric considerations. It is T = - + (dS+dR) cosi V2 Vi (L Lcrit). (2.5) Β£ is the source-receiver range measured parallel to the refracting interface and d and dR are perpendicular distances from S and R to this interface. The critical refraction angle i is given by the usual expression: sin i = vi/v2. Of course, the head wave exists only for ranges exceeding the critical range given by Lcrit = (d + dR) tan i. R (2.6) Fig. 2.3. Head wave raypath. The plane of this diagram is perpendicular to the refracting interface and thus is not necessarily vertical. S L ds dR P Q 14 The head wave traveltime formula is more useful to the geophysicist if it is expressed in terms of the horizontal offset X. A simple analytical derivation of the desired relation is given in this section. An advantage of this method is that it does not require visualization of the three dimensional geometry of the problem. Starting with the general equation for the dipping plane, it is easy to demonstrate that the perpendicular distances from S and R to this plane are given by d = d β rs n and dR d β rj . n, respectively. Hence dR = ds+(rsβrR).n. Substituting from equations (2.2), (2.3), and (2.4) and reducing yield = dsβXsinq5cos(aβ6). Now define an angle 6 (βTr/2 < S < ir/2) as follows: sinS sin q5 cos(a β 8). Then, the above expression becomes dR = dsβX sinS. (2.7) It is stressed that the angle S is not the apparent dip of the refracting interface along the profile line azimuth a. Apparent dip -y is related to true dip Γ§tβ via the following equation: tan = tanq5 cos(aβO). (2.8) Obviously, S βy. The difference between S and y arises from the fact that the former is measured in the raypath plane, whereas apparent dip is measured in a vertical plane. For small values of true dip, S -y. The range Β£ is equal to the distance FQ in Figure 2.3. Hence 15 L2 = = (rS+dSn)β(rR+dRn) 2 = (rS_rR)+(dS_dR)n2 = IrsβrI?I2+2(dSβdR)(rSβrR)n+(dsβdR) = = x2 β x2 s. Thus L = X cos8. (2.9) Equations (2.7) and (2.9) are the desired expressions. It is evident that these formulae could have been derived by purely geometrical reasoning based on the raypath diagram of Figure 2.3. However, the relationship of the angle S to the proffle line azimuth and the orientation angles of the refracting plane would not have been easy to ascertain. Substituting expressions (2.7) and (2.9) into equation (2.5) and reducing yield an expression for the head wave traveltime: T = sin(icS) + 2dscosi (XXcrit). (2.10) An expression for the critical offset distance X is obtained by similiar manipulations. Evaluating equations (2.7) and (2.9) at the critical distance and then substituting into (2.6) gives slnic / = 2d, . . 2.11 cos(z β 6) The crossover distance is defined as the particular source-receiver offset at which direct wave and head wave arrival times are identical. Setting T in (2.10) equal to X/vi and solving for X yields cos z X035 = 2d . . (2.12)1 β sin(z β 8) Fig. 2.4. Normalized apparent refractor velocity Va/V2. Each curve refers to a specific value of the interface dip angle . The P-wave velocity ratio is vl/v2 = 3/5. 16 Finally, the traveltime formula (2.10) indicates that refracted arrivals propagate along the receiver spread with an apparent velocity Va given by V1 Sifl 2 Va . β’ V2. (2.13) sin(i β 8) sin(z β 8) The variation of apparent velocity with profile azimuth is depicted in Figure 2.4 for the specific case where vl/v2 = 3/5. Equations (2.10) through (2.13) indicate that the three dimensional refraction formulae are straightforward generalizations of those appropriate for the 2D problem: the angle S replaces the true clip S in the relevant expressions. It is emphasized again that S is not the apparent dip of the subsurface interface along the profile line. However, for a gently dipping refractor, the practical difference between the two is small. 3.0 2.5 > >t5 1.0 0.5 0 100 200 300 aβS (deg) 17 2.4 Traveltime inversion Traveltimes recorded by a set of refraction proffles can be inverted to recover the 3D attitude, true velocity, and depth of the refractor. Let the measured slope and intercept time of the th such traveltime curve be m and r, respectively. From equation (2.10), these are related to the assumed earth model parameters via sin(i β S) 2d. cos m = = , (2.14a,b) vi Vi with: sinSj = sinqScos(ajβ9), (2.14c) ds d β sin q (zs cos 8 + ys sin 8). (2.14d) (zs, y.g) are the position coordinates of source i and a: is the azimuth angle of proffle line i. The perpendicular depth to the interface below source i is designated How many refraction lines are required in order to successfully invert for the earth model parameters? Assuming that the overburden velocity vi is known (perhaps from borehole data or traveltime analysis of the direct arrivals), then the earth model is defined by the four parameters (v2, q, 9, d). Intuition suggests that two refraction profiles will yield the four data (two slopes and two intercepts) required to solve the problem unambiguously. Indeed, if the line index i is set equal to 1 and then 2, the system (2.14) above becomes a set of 8 equations in 8 unknowns. Since the equations are nonlinear in the unknowns, a definitive statement on the existence and uniqueness of a solution cannot be made. However, it can be demonstrated by algebraic techniques that, in many cases, two refraction proffles are sufficient to solve the problem. The method proposed by Russell et al. (1982) requires three lines to obtain the re fractor velocity, attitude, and depth. They utilize the slopes of the three traveltime curves (mi, m2, m3) and one intercept time (ri) as measured data in the inversion. For the case where the index i runs from 1 to 3, system (2.14) represents 12 equations in 10 unknowns; 18 hence, it would seem that redundant information exists. Although redundant data are al ways valuable in any practical inversion scheme, the theoretical minimum set of conditions under which the inversion is possible are of interest here. In the following analysis, only two refraction profiles are employed. First, define the angles P1 and 112: sinβ(mlvl), sinβ(m2vl). (2.15a,b) Next, the position of source 2 is expressed in terms of the coordinates of source 1: Xs2 = xs, + H cos/3, ys2 = ys1 + H sin /3. (2.16a,b) /3 (0 /3 < 2ir) is the azimuth of source 2 relative to source 1 (measured positive in the clockwise sense from geographic north) and H is the horizontal distance between the two source locations (see Figure 2.2b). Using these expressions, the angles 81, 8 and the depths ci, ds,, cl2 can be quickly eliminated from the system (2.14). A reduced system is obtained consisting of three nonlinear equations in the three unknown angles i, 4, and 0: sin(i β 111) = sin 4 cos(ai β 8), (2.17a) sin(i β 112) = sin 4 cos(a2 β 8), (2.17b) vl(T1 β T2) = 2H cos i sin cos(/3 β 0). (2.17c) 24.1 Special cases Three particular cases of the above system are examined before a more general solution technique is described in the next section. 1) Coincident sources. If the two sources occupy the same position (H = 0), then equation (2.17c) above reduces to 0 = 0 (coincident source locations imply identical intercept times). 19 The two remaining equations contain three unknowns and cannot have a unique solution. Hence, two refraction proffles emanating from the same source location supply insufficient independent information for a successful inversion. The split recording spread is a particular example of coincident-source refraction profiles (Johnson, 1976). 2) Parallel profiles. If the azimuths of both refraction lines are the same (cr1 = cr2), then equations (2.17a) and (2.17b) above are not independent (identical azimuths imply identical slopes of the traveltime curves and thus pi = pa). Hence, the parallel profile recording configuration is not adequate to solve the problem either. 3) Anti-parallel profiles. Suppose that line 2 is recorded in a direction opposite to that of line 1 (cr2 = cr1 Β± jr). Then equations (2.17a) and (2.17b) can be solved immediately for the critical refraction angle: = 1U1 + P2 (2.18) The refractor velocity is obtained via V2 = v1/ sin c. This procedure for calculating V2 is identical to that used with the classical reversed spread in the two dimensional situation. It is clear that it is generally valid for anti-parallel proffles (colinear or otherwise) recorded over a three dimensional dipping interface. With the critical angle i determined, the system (2.17) can be solved for the remaining unknowns. The azimuth of line 1 is designated by a. Straightforward, but tedious, algebraic manipulation then yields tanO β β (r β7-2) cosa + H(mi βm2) cos/3 (219) β (Ti β T2) sin a + H(mi β m2) sin/3 sin 4β = F J(i. β 7-2)2 + 2(ri β T2) H (ml β m2) cos(a β /3) + H2(mi β m2)2, (2.20) where the multiplicative factor F is given by F . . (2.21) 2H cos i sin(a β /3) 20 A minor ambiguity associated with selecting the proper branch of the inverse tangent func tion in equation (2.19) is easily resolved by ensuring that the angle 8 satisfies the original system (2.17). After 0 and Γ§b are determined, one of the equation pairs (2.14b,d) can be solved to yield the distance d from the reference point 0 to theinterface. The solution for the two orientation angles given by equations (2.19) and (2.20) contains an important indeterminacy when the azimuth angle /3 of source 2 equals a or a Β± ir. In this case, the recording geometry consists of two colinear, anti-parallel proffies. The identity of reciprocal times recorded at the two inline shot positions (i.e., mlH+Ti m2H-Hr)implies that each equation reduces to 0/0. Hence, two colinear anti-parallel proffles are inadequate for defining three dimensional refractor attitude. The reversed spread is a typical example of this configuration. As long as there exists some perpendicular offset distance between the two lines recorded in opposite directions, the inverse problem is well posed, at least in a theoretical sense. The latter data acquisition geometry is a fairly common arrangement for many land and marine seismic surveys. Arbitrary line azimuths A combination of algebraic and algorithmic techniques yields a solution of the system (2.17) for the case of arbitrary proffle line azimuths al and a. The special situations treated in the previous section are excluded. Eliminating 8 and Γ§S results in a quadratic equation in tan i: A tan2i+ B tani + C = 0. (2.22) The coefficients A, B, and C depend upon measured quantities and recording geometry as follows: A V1T1T2 f = 2H sina β B = β cos sin(a2 β 3) + cos 1β2 sin(ai β /3), (2.23b) C sin sin(a2 β /3) β sin2 sin(ai β /3) + A. (2.23c) 21 In the case of anti-parallel proffles the constant A vanishes and tan i is obtained by solving a linear equation. It is easy to verify that the expression for the critical angle i given by the previous equation (2.18) is reproduced. In general, solution of equation (2.22) is via the quadratic formula: tan βB + sgn(B) β 4AG (2.24) 2A where the sign of the radical has been chosen to yield an indeterminacy for A = 0. Ap plication of LβHopitalβs rule then indicates that (2.18) is recovered as A β* 0. After i is determined, simply back substitute sequentially into the expressions for the remaining unknowns. This is done numerically rather than symbolically since the formulae rapidly become unwieldly. Hence tan9 = β H cosi sin(ic β p.i) cos/3 β [vl(T1 β r2)/2] C05a1 (225) H cosi sin(i β iβ) sin3 β [vl(r1 β βr2)/2] sinai An analogous (and completely equivalent) formula involving 2 and a2 could be used. Again, care must be exercised to ensure that the angle 8 obtained by inverting (2.25) satisfies the original system of equations. Any one of the expressions (2.17a,b,c) can then be solved for the true dip angle q5. Finally, depth d is determined as previously described. 2.5 Traveltime isochrons Heretofore, the analysis is concerned with an inversion method for critically refracted traveltime data observed on two line profiles. An alternate approach to interpreting three dimensional refraction data entails recording arrival times on an areal grid of receivers. The measured data are posted and contoured on a map display of the receiver array. Charac teristic patterns observed in the contoured data are used to qualitatively infer subsurface structural style. If the traveltime data are sufficiently precise, a quantitative estimate of various structural or stratigraphic parameters (dip and strike angles, locations of faults or 22 edges, magnitude and sense of fault throw, refractor velocity, etc.) is possible (Dunkin and Levin, 1971; Palmer, 1986, p. 246-249). For a fixed traveltime T, the head wave wavefront occupies a three dimensional spatial locus consisting of a frustum of a cone. The intersection of this locus with the recording plane z = 0 is defined to be a traveltime isochron. A closed form expression for a head wave traveltime isochron in polar coordinates is easily obtained by solving equation (2.10) for offset X as a function of the azimuthal angle a: vi(T β r) X(a) = . . , (2.26) sin{i β 6(a)] where T is the intercept time (r = 2d cos i/vi) and 8(a) sin1 [sin Γ§ cos(cr β 9)]. Straightforward algebraic manipulation reduces this expression to a quadratic form in the distance X: X2 1β (Sin )22 β9) β 2X v1 (T β r) COS cos(a β9) β v (T β r) 2 = slnic sin i sinic (2.27) The mathematical structure of this expression is identical to that of equation (A2) in Ap pendix A. Hence, a locus of fixed head wave arrival time T is a conic section (either an ellipse or a hyperbola) with principal axes parallel to the dip and strike directions of the refracting interface. Furthermore, the parameters defining the geometry of the curve can be determined by equating coefficients in (2.27) and (A5). Solving for the eccentricity e, semimajor axis a, and the origin-to-center distance l yields β sin q β vi(T β r) sin i cos β vi(T β r) cos i sin qS e β , a β . 2 .2 β β’ 2 .2 sin z sin z β sin sin z β sin (2.28a, b, c) These three expressions are consistent with equivalent formulae derived in an entirely dif ferent manner by Dix (1935). Two distinct cases must now be examined. 23 Case 1: > q1. If the critical angle exceeds the interface dip angle, then the eccentricity of the conic section is less than one (e < 1). Hence, traveltime isochrons are effipses that encompass the source position. However, the shotpoint is not located at a center or a focus of any effipse. Rather, the center is displaced in the updip direction by the distance I, and the nearest focus is displaced downdip by the amount a e β 4. For a horizontal interface, the above expressions reduce to e = 0, a = v2(T β r) and l = 0. In this simple situation, head wave traveltime isochrons are concentric circles (centered at the source position) with radii that increase linearly with T β T. Case 2: i < . In this situation, the eccentricity exceeds one (e > 1) and isochrons are hyperbolae. This case may arise where there is a strong velocity contrast between overburden and refractor, together with a large dip on the intervening interface. Head wave traveltimes actually decrease with offset distance in a wedge of azimuths centered about the updip direction (i.e., T β r < 0 for c bounded by Β± cos1(1/e)). Elliptic and hyperbolic isochrons for head wave traveltime are illustrated in Figure 2.5. These contours are associated with actual head wave arrivals only if the source-receiver offset exceeds the critical distance given by equation (2.11). The curve defined by Xt(a) delimits the precritical offset zone, within which no head waves are observed. Expression (2.11) easily reduces to [i β (sin 1) 2 cos(cz β 8)] β 2Xcrit [β2ds tan2 i sin 4β cos(cr β 8)] β [2d tan c] 2 = cos z (2.29) Hence, the limit of the precritical offset zone is also a conic section. By correspondence with equation (A4), the parameters defining this conic are sin 4β 2d sin i, cos cos 4β β2d5 sin2 i sin 4β e = , a = 2 β’ 2 , = 2 β’ 2 cos i cos z β sin 4β cos i β sin 4β 24 400 200 E -c cβ20 β400 400 200 0 -c 0 c β200 β400 β400 400 easting Fig. 2.5. Head wave traveltime isochrons. In each panel, the isochrons (contour interval = 25 ms) are represented by solid curves, the critical offset distance by a dotted curve, and the crossover distance by a dashed curve. The source (asterisk) is located at the origin, and the straight line indicates the outcrop locus of the refractor. (a) Earth model: vi = 1500 m/s, v2 3000 m/s, q. = 15Β°, 0 = 45Β°, l = 100 m. Isochron eccentricity = 0.518. (b) Earth model: v = 1500 m/s, v2 = 5796 m/s, = 30Β°, 0 = 45Β°, h = 100 m. Isochron eccentricity = 1.932. 0 (m) 25 400 200 E 0 -c - β200 β400 400 (m) Fig. 2.5. (c) Earth model: v1 = 1500 m/s, v2 = 1655 m/s, g = 300, 6 = 450, h = 100 m. Isochron eccentricity = 0.552. If < ir/2 β i, the eccentricity is less than one and the curve defined by (2.11) is an ellipse surrounding the source; its center is shifted by an amount l in the downdlip direction. However, if the interface dip angle exceeds the complement of the critical angle (4> r/2βi), then the critical offset curve is a hyperbola. This occurs, for example, in a high dip, low velocity contrast environment where rays critically refracted downdip never return to the surface. The limit of the precritical offset zone is displayed in Figure 2.5 as a dotted curve. Fictitious arrival times within this curve are obtained by extrapolating or phantoming real traveltimes. Figure 2.5a illustrates elliptic traveltime isochrons ( > 4) and an elliptic critical distance curve (i + 4 < ir/2). In Figure 2.5b, the isochrons are hyperbolic (i < 4) and the critical offset is still elliptic. Finally, in Figure 2.5c, the traveltime contours are β400 0 easting 26 elliptic (ia> q) and the critical distance curve is hyperbolic (i + q> Tr/2). The remaining case (hyperbolic isochrons and hyperbolic critical offset) is not illustrated. The crossover distanceX03(a), computed from equation (2.12), is depicted by a dashed curve in each panel. Since expression (2.12) cannot be put into the form of equation (A4), this curve is not a conic section. 2.6 Inversion of traveltime isochrons If the geometric parameters of a traveltime isochron can be determined from head wave arrival times observed on an areal recording array, then it is possible to invert for the earth model parameters. The principle axes of an effiptic or hyperbolic isocliron are parallel to the dip and strike directions of the subsurface refracting horizon. Hence, the azimuth angle & can be found. Furthermore, solving equations (2.28) for sini, sin, and vijT β r in terms of e, a, and yield 1 1a2e2 β 12 /a2e2 β sini = β 2 β sin = 2 (2.30a,b) vilTβri = /(i βe2)(a βl). (2.30c) For an ellipse, these expressions are all well defined because e < 1, a > l, and ae > l. However, the relations remain valid for a hyperbola where e> 1, a < l1, and ae < lI; the argument of each square root is still positive. Equations (2.30a,b) indicate that the critical and dip angles can be determined from the geometric parameters of an observed traveltime isochron. If the overburden velocity vl is known from ancillary data, then(2.30c) can be solved for the intercept time r. The normal depth to the refractor below the shotpoint is then calculated via d = v1T/2 cos i. However, an alternative approach considers the intercept time to be a measured quantity. The basic traveltime equation (2.10) indicates that intercept time is invariant with respect to recording profile azimuth c. Thus, a value for T can be obtained by linearly extrapolating 27 to zero offset the arrival times observed along any radial profile emanating from the source. Equation (2.30c) is then a single equation in the unknown layer velocity v1. Solving for v1 gives ___ ___ ___ ___ ___ ___ βe2)(a βl) Vj = eITβrI . (2.31) Hence, the overburden velocity can be calculated from the refraction data alone! With v1 determined, the refractor velocity is obtained via v2 = vi! sin i. Interestingly, the above solution method requires the critically refracting interface to have nonzero dip. Otherwise, the eccentricity e and origin-to-center distance 4 vanish, and formula (2.30c) becomes inde terminate. Overburden velocity v must then be obtained from independent measurements. Practical issues related to fitting a conic section to a finite set of spatially sampled arrival times are not addressed here in any detail. Rather, the intention is to demonstrate the existence of a mathematical solution to the specified inverse problem. In principle, the situation is simffiar to fitting a straight line segment to a set of error contaminated traveltimes observed along a line profile. However, there are obvious practical differences. Since a conic curve is nonlinear in its defining parameters, an iterative fitting procedure may be necessary. Finally, it may be possible to fit several conics with different isochron times T to the measured data. Statistically robust estimates of the earth model parameters can then be obtained by using some form of averaging. 2.7 Conclusion Although the earth model examined in this chapter is quite simple, there are undoubt edly many situations where it constitutes an adequate representation of geophysical reality. Also, the model is very useful as a pedagogical tool. Analysis of head wave propagation within the model yields closed form mathematical expressions for traveltime, critical and crossover distances, intercept time, and apparent refractor velocity. These expressions are generalizations of more familiar formulae that apply in one and two dimensional situations, and are valuable for forward modeling purposes. 28 Algebraic solutions of the traveltime inverse problem for the same earth model are also derived. There are numerous practical recording geometries where the analysis of critically refracted arrival times recorded on only two unreversed spreads can yield the three dimen sionai attitude, true velocity, and depth of a plane subsurface horizon. The well known method for determining refractor velocity in a 2D situation (measurement of apparent ve locities on colinear forward and reverse profiles) applies directly to the more complicated 3D model. However, it is not necessary that the proffles be colinear; anti-parallel profiles are sufficient. This knowledge can be quite useful for the interpretation of apparent velocity - data obtained from a 3D dipping interface. Extension of these methods to include multilay ered earth models will be valuable for defining the minimum requirements for a successful inversion of three dimensional refraction traveltime data. Finally, study of head wave arrivals recorded in an unconventional data acquisition geom etry (i.e., an areal receiver array) reveals an interesting possibility for estimating overburden velocity from the refraction data alone. At present, this is primarily a theoretical result, but it does provide some justification for the numerical inversion results obtained in Chapter 4. 29 CHAPTER 3 HEAD WAVE TRAVELTIMES IN A THREE DIMENSIONAL MULTILAYERED EARTH 3.1 Introduction There is a notable paucity of papers on the subject of head wave propagation in three dimensional, multilayered earth models. Chander (1977b) examines a model consisting of uniform velocity layers separated by plane interfaces with arbitrary strike and dip, and de scribes a method for calculating head wave traveltimes between specified source and receiver positions on a horizontal surface. If an array of receivers is colinear with the source, the head wave arrival time curve is a straight line. Hence, if two points on this line are established, then arrival times at all offsets can be determined simply by drawing the connecting straight line. Chander (1977b) locates the two initial points via raytracing techniques. Chanderβs work is purely numerical and does not provide much insight into the depen dence of head wave traveltime on the parameters that define the earth model. Moreover, it is restricted to conventional data acquisition geometries. Buried sources and/or receivers as well as nonproffle recording geometries require a more general treatment. Dieboldβs (1987) recent work constitutes the seminal contribution on this topic. He also considers a three dimensional multilayered earth, but derives traveltime formulae for reflected and critically refracted waves. These formulae are logical extensions of familiar traveltime expressions that are appropriate for one and two dimensional layered models. Thus, they offer the possibility for extending several known traveltime inversion techniques to accommodate 3D planar structure. Unfortunately, Dieboldβs derivations are very ambiguous. Furthermore, his generalization to arbitrary source-receiver geometries yields an incorrect traveltime formula. Finally, he does not present a numerical technique for actually computing the traveltimes. These deficiencies are addressed in this chapter. Nevertheless, Diebold (1987) should be credited with an original contribution to traveltime analysis for this particular class of earth models. 30 This chapter provides a rigorous derivation of the three dimensional head wave traveltime formula. A related expression for the traveltime of a reflected wave propagating in the same earth model is obtained as a byproduct of the analysis. The mathematical proofs of the formulae are simplified by using a novel form of Snellβs law of refraction and reflection. Various generalizations of the basic traveltirne equation extend its applicability to arbitrary 3D recording geometries and/or mode-converted waves. Finally, a rapid numerical method for computing the arrival times of critical refractions is presented, and is illustrated with simulated examples from shallow refraction exploration and vertical seismic profiling (VSP). 3.2 Earth model Consider an earth model consisting of a set of homogeneous and isotropic layers bounded by plane interfaces. In general, each interface may possess a three dimensional dipping attitude. The interface of the model is illustrated in Figure 3.1. 0 is the origin of. a right-handed, rectangular Cartesian coordinate system with orthonormal basis triad ijk. The xy plane is defined to be the horizontal plane and the depth coordinate z increases in the downward direction. The locus of plane interface i satisfies the equation rβ’n = dj, (3.1) where n is a unit vector normal to the interface and dj is the perpendicular distance from 0 to the interface. Figure 3.1 indicates that n is conveniently described by two interface orientation angles: n; = (sin Γ§fj cos 8)i + (sin Γ§t.j sin 9)j + (cos g5)k. Γ§j (0 <q <ir/2) is the dip angle and O (0 <O <2r) is the azimuth angle of the interface. If the +x and +y axes are taken to point toward geographic north and east, respectively, then the interface strike angle is 8 + r/2 (modulo 2ir). Although these angular coordinates 31 y 0 Fig. 3.1. Orientation of the iβ interface of a multilayered earth model. are descriptive, a certain compactness in notation is achieved by specifying n2 in terms of its Cartesian components: n = nj,i + ni,,i + x ni z vi +1 with n + m, + = 1. This convention is followed in the sequel. 32 Solving equation (3.1) for z as a function of z and y yields the vertical depth of interface zQc, y) z(O, 0) β β y [3], (3.2) Th:,z where z(O, 0) d/nj, is the vertical depth of the th interface below the coordinate origin. The surface (not necessarily horizontal) is interface 1, and subsequent interfaces are numbered sequentially in the downward direction. Interface i overlies layer i. The vertical thickness of layer i is defined to be h(z,y) zj+1(x,y) β zj(x,y). Thus h(x,y) = h(0,O) + [Th β hh+ij β (3.3) βi,z i+1,z fl,z i+1,z where h(0, 0) z+i(O, 0) β z(O, 0) is the vertical thickness of the layer beneath the coordinate origin 0. Finally, the seismic wave propagation speed assigned to layer i is given by v. This may be either the compressional wave speed cx or the shear wave speed f3. This fiexibffity allows the resulting traveltime equations to apply either to P, S, or mode-converted waves. 3.3 Raypath geometry Initially, the analysis is restricted to the case where both the source and receiver are located on the surface. Generalization to an arbitrary data acquisition geometry is straight forward and is given in a later section. The horizontal coordinates of the point source S and point receiver R are (as, Ys) and (R vi?), respectively. Their vertical coordinates are easily obtained from equation (3.2): z = zl(zs,ys) and ZR = zl(XR,YR). In order to facffitate computation of the traveltime, the total head wave raypath is divided into three major portions: the downgoing, critically refracted, and upgoing paths. In Figure 3.2, these correspond to raypath segments SF, PQ, and QR, respectively. The propagation time along each portion is calculated, and then all three are summed to obtain the surface-to-surface head wave traveltime. 33 Fig. 3.2. Schematic representation of the raypath of a wave critically refracted on interface k of a multilayered earth model. S and R denote a surface source and receiver, respectively. In the 3D situation, the raypath is not confined to a plane. Within each layer, the raypath is a straight line segment. On the downward portion of the raypath, the propagation direction within layer i of the wave critically refracted at subsurface interface k is described by the unit vector pj: Pik = Pik,z + PiJc,.i + S R . . . . P 0 Similarly, the upward propagation direction within layer i of the wave critically refracted 34 from interface k is specified by another unit vector qj: qk = q,zi + qtk,ij + q:,zk. The complete head wave raypath is described by the set of unit vectors Pik and q (i = 1,2, . . . , k β 1) together with a critically refracted propagation direction pjj, = qkk At interface i in the overburden, the wave is refracted in accordance with Snellβs law. The situation for the downgoing wave as it encounters the iββ interface from above is depicted in Figure 3.3a. The plane of this diagram is the plane of incidence defined by the incident propagation direction Piβ1,k and the interface normal n. Snellβs law of refraction consists of the following two conditions: (i) the unit propagation vector of the transmitted ray (pjk) is contained in the plane of incidence, (ii) sin Pjvj_1 = sin v/v, where and v are positive acute angles measured from the in terface normal to the incident and transmitted propagation directions Piβ1,k and Pik, respectively. Both conditions are contained in the single vector equation fli )< Piβi,k = n X Pik (3 4) Vj_1 V The vector formed by these cross products points out of the plane of the diagram in Figure 3.3a. In component form, equation (3.4) is 35 Fig. 3.3. Noncritical refraction of the raypath at interface i in the overburden. (a) Down- going raypath. (b) Upgoing raypath. The unit vector n is normal to the interface. v-1 Pi-1,k a Interface I ni P1k iflterfa I qi-1,k b vββ1 qik ni 36 1 r 1 ir 1 {fli,Pβ1,k,z β flizPi_i,k,yj β [fli,yPik,z β fli,zPik,j (3.5a) v_1 Vi [ni,zPi_i,k,z β fli,zPi_1,k,z] = [rti,zPilc,z β fli,zPik,zj, (3.5b) V2_ Vt ii 1 ir 1 [fli,zPiβ1,k, β fli,vPiβ1,k,zj β [fli,zPik,i, β flz,Pik,zj. (3.5c) Vj_1 Similarly, when the upgoing wave encounters interface i from below, Snellβs law in the form flj X qβ1,k ni x (3 6) Vj_1 V holds (see Figure 3.3b). In this case, the angles 1ti and v exceed ir/2 radians. The component form of expression (3.6) is analogous to equations (3.5). The three dimensional statement of Snellβs law of refraction given by the above expres sions is quite different from the form typically used in raytracing applications (Sorrells et al., 1971; Shah, 1973; Chander, 1977b). However, it can be demonstrated that these expressions are equivalent to the raytracing formulae (equations (3.27) and (3.28) below). The value of the current formulation is that it leads to a substantial simplification in the mathematical proof of the traveltime equations. 3.4 Traveltime derivation 3.4. 1 Downgoing traveltime An expression for the traveltime increment of the downgoing wave as it traverses layer i is derived first. Let the position vectors r: and r:+I denote the intersection points of the downgoing ray with interfaces i and i + 1, respectively. Then rj+1 = r + 13 Pik, where 37 4 is the length of the (straight line) raypath segment within layer i. Solving for 4 gives 4 Pk (r+I β rj. In component form, this expression is 4 = pjk,z(Xi+1 β x) + pik,y(yi+1 β y) + Pjk,z(Zi+1 β z). The traveltime increment is obtained by dividing the path length segment 4 by the layer velocity v: ii Pik,z Pik,y Pik,z ti = β = ββ(xi+i β x) + β(yi+1 β yi) + β(zj+i β z). vi v vi vi The vertical (z) coordinates of the intersection points can be expressed in terms of the horizontal (x, y) coordinates by using the equation for a dipping plane interface. Equation (3.2) yields rni,zl Fβi,yl z zj(xj,yj) = zj(O,O)βajiβ-βi βyiββ--, LThi,zJ Lfli,zJ β rni+1,zl rni+i,y Zj+1 = zj+1(xj+1,yj+l) = zi+1(OO)βxi+1[ ] Yi+4 Hence 1i,1 f2i,v1 1βi+l,z1 1βi+i,y1 β = h(O, 0) + z βj + β x+1 [ j β Yi+1 [Th:,z where h(0, 0) = z+i(0, 0)β z(0, 0) has been used. Substituting this result into the equation for t and grouping terms yields the required expression for the traveltime increment: = h(O,O)pk, + xi+l {fli+1zPikz fli+1zPikz] + Yi+1 {fi+1zPik hui+1Pikz] Vi Thi+1,z V nβi+l,z Vi β [,z1 β fli,Pik,z1 β L [,zPik,y i,yPik,z vi -I L v 38 The total downgoing traveltime is obtained by summing all of the individual traveltime increments t, i = 1,2,. . . , k β 1: kβi ___ ___ ___ __ Idown = Vi + 2i+i [fli+izPikz β7Zβi+i,zPik,z] + Yi+i [fl-i+izPik β Thi+iPikz] i1 βj+1,z β i [flizPikz β i,zPik,z j + [flizPikv β fliPikz] i=I (Thi,z V2 V fli,z V V The sum involving aji and Yi-i-1 is now reindexed and combined with the other sum. The result is kβi β Γ§ βdown β vi + Xk(flk,zPki,k,z β Thk,zPkβi,k,z) + Yicfrβic,zPk_i,ic,7,β Thk,11Pkβi, ,z) Vk_iflk,z β Xi(Th1,zPik, β fli,zPik,z) + y1(fli,zPiic,y β fli,yPlk,z) Vimβ,2 + i_ {TiizPi_ikz β i,zPiβi,k,zj β [flizPiJcsz β i,Pik,z] . ni,z I Vi_i Vi :=2 + {flizPi__1k i,yPiβi,k,z] β {flizPik β i,yPik,z] i=2 Vj_1 vi 39 At interfaces 2,3,. . . , kβi in the overburden, Snellβs law of refraction must apply. Equations (3.5a,b) then imply that all terms in the summations involving x and Y: vanish! The downgoing traveltime from (x1,yl) = (xs,ys) to (zk,yk) (xp,yp) reduces to Icβi ___ ___ ___ __ .Ldown vi + XP(Thk,zpk_1,lc,z β flk,zPkβ1,k,z) + yP(Thic,zpIcβi,k,y β Thk,yPlcβi,k,z) Vk_iflkz β ZS(fli,zpik,z β mi,zpilc,z) + YS(fli,zplk,y β fl1,yplk,z) (37) Vlfll)Z 3.4.2 Upgoing traveltime The traveltime along the upward propagating portion of the total raypath is derived by similar techniques. Snellβs law in the form (3.6) is used at each interface in the overburden. The result is kβi βΓ§β ___ βupβ Ld i=i vi β XQ(flk,zqk.i,k,z β flk,zqk_1,k,z) + yQ(nk,zqk_i,k,y β nk,yqk1,k,z Vk.lflk,z + ZR(Th1,zqik, β ni,zqlk,z) + yR(ni,zqlk,y β n1,qik,z) (38) vβnβ,z 40 3..3 Critically refracted traveltime The final traveltime increment needed for the derivation corresponds to the critically refracted segment of the total raypath. Let the position vectors rp and rq refer to the intersection points of the downgoing and upgoing portions of the raypath with interface k, respectively. The critically refracted raypath segment is a straight line connecting these two points (and thus lying entirely within the plane of interface k). Then rq = rp + tk Pick and hence l = Pick . (r β rp). The propagation time along this path length is Plck,z Pkk,y Pkk,z tic = β = (x β zp) + (y β yp) + (z β zp). Since points P and Q reside on the same plane interface, equation (3.2) yields rk,z1 ZQZp = Zk(Q,yQ)-Zk(Zp,yp) = β(q---a)-ββj fk,z k,z The expression for the critically refracted traveltime increment then reduces to β β (xq β Xp)(rLlc,zpkk,z β k,zPkk,z) + (Y YP)(k,zPkk,y Thk,yPkk,z) tic β Vkflk,z (3.9) 3.4.4 Total traveltime The total surface-to-surface traveltime of the wave critically refracted on interface k is obtained by adding the traveltime contributions along the downgoing, critically refracted, and upgoing raypath portions: Ttotai = Td0 + + Summing expressions (3.7), (3.8), and (3.9) yields 41 kβi Tk(XS,YS,aR,YR) β qiie,) vi i=i + [xi(n1,zqβ1z β n1,q1k,z) + yR(nl,zqik,y β m1,q1k,z)] Vifli,z [ ZS(flI,zPik,z β fli,zPik,z) + YS(nl,zpik,y β fliYP1kz)] vlnl,z + F(XP,YP,XQ,YQ). The quantity F depends on the horizontal coordinates of the two points of critical refraction and is given by F(xp,yp,XQ,yQ) = xpF 1 1 I (k,zPkβ1,k,z β Thk,zpkβ1,k,z) β (βk zPkk,z β ThkzPkkz)] k,z Lvk_1 Vk Fi 1 + I (flk,zPkβ1,k,y β flk,yPkβ1,k,z) β (Thk,zPkk,y β flkPkkz)] k,z [vk_1 Vk _ 1 I (n,zq_1,,z β nk,zqk_1,k,z) β β ThkzPkkz)] k,z [Vk_i Vk 1I (βk,zqkβ1,k, β nk,yqk_1,k,z) β (k,zPkk,y β flkiPkkz)]. k,z [vk_I Since the wave is critically refracted at interface k, the propagation direction vectors Pkk and kk are identicai. Then, requiring Snellβs law to be satisfied at points P and Q results 42 in F(xp, yp, XQ, y) = 0. The final formula for the surface-to-surface traveltime of the head wave becomes kβi β h1(O,O)(p,, β qik,z) Tk(xs,ys,XR,yR) = i=1 $ + xR(n1,zqik, β nI,zqlk,z) + yR(n.1,zqik,y fll,yqik,z) Vifli,z β ni,zplk,z) + YS(nβi,zpik,y β Thi,yPiIc,z) (3.10) [ vini, 34.5 Variants of the basic formula Obviously, if the source or receiver is located at the coordinate origin, then the traveltime formula (3.10) simplifies considerably. Another simplification arises with a. horizontal surface (ni, = = 0, n = 1). Expression (3.10) reduces to Tk(xs,Ys,XR,YR) = h(o,o)(Pik β qik,z) + + yRqik,y)β (ZSP1k,z + YSP1k,y)] (3.11) This is equivalent to Dieboldβs (1987) surface-to-surface traveltime formula. Conventionally, seismic refraction traveltime is presented as a function of the source- receiver offset distance. The current expression is easily converted to this form by specifying the receiver position in terms of an offset distance X (X 0) and an azimuth angle P (0 < Eβ <27r) relative to the source. The receiver coordinates are XR x+XcosβIβ, YR = ys+Xsin. 43 Substituting these expressions into equation (3.10) yields kβi I V T Γ§β h(0, 0)(Pk,z β qik,z) 1kkXs,yS,,) = Thi,z(pik,x β qik,z) β 21,z(P1k,z β qlk,z) β 2s ViThi,z fll,z(pik,y β qlk1y) β Th1,y(Pik,z β qik,z) β ys Vifl.1,z + cos (n1,zq1, β Th1,zqlk,z) + Sfl (ii,zqiic,y β ni,yqik,z) (3.12) Vlfli,z Note that X is the horizontal distance between source and receiver; the actual distance is larger since it is measured within the plane of interface 1. It is straightforward to demonstrate that the true source-receiver distance is L XsJl + tan2 li cos2(β 8), where 4βi and 6 are the dip and azimuth angles of the surface. Equation (3i2) is an extension of the common βslope and interceptβ traveltime formula for head waves to three dimensional, multilayered earth models. For the particular case of a model with only two layers and a horizontal surface, it can be shown that this equation reduces to the closed form expression (2.10) derived in the previous chapter. This serves as an important check on the correctness of the general formula. However, the proof entails some cumbersome algebra, and hence is relegated to Appendix B. Finally, consider the specialization of the general traveltime formula to a two dimensional earth model. In the model parameterization utilized here, all of the y components of the interface normal vectors vanish (equivalently, all interface azimuth angles are restricted to = 0 or = ir). Additionally, the recording profile must be oriented perpendicular to the strike directions of the subsurface horizons. Hence, I1 = 0 or r also. If this second condition 44 is not satisfied, the wave propagation vectors p,j and qjj are not confined to the xz plane and a full three dimensional treatment is necessary. Since many previous investigators have assumed a horizontal surface, equation (3.11) is used as the point of departure for the analysis. Setting s = YR = 0 gives Tk(XS, XR) = hi(O, O)(P:lc:z β qik,z) + qik,z XS Plk,zj (3.13) - where +Pkz = + qk = 1. This expression is compatible with analogous equations developed by Diebold and Stoffa (1981) and Diebold (1987). If the source is located at the coordinate origin, then (3.13) is also consistent with earlier two dimensional head wave traveltime formulae published by Dooley (1952), Adachi (1954), Ocola (1972), and Johnson (1976). All of these investigators prescribe vertical layer thicknesses, either beneath the origin or the shotpoint. In contrast, Ewing et al. (1939) and Mota (1954) measure thickness normal to the basal interface bounding a layer. Hence, their traveltime equations, although designed to treat an equivalent situation, differ in mathematical detail. 3.5 Generalizations of the traveltime formula Heretofore, both the source and the receiver have been restricted to the surface. More versatile formulae are needed to model data acquisition geometries with buried sources and/or receivers. These situations arise in surface-to-borehole, borehole-to-surface, and borehole-to-borehole seismic experiments, as well as with placement of sources and/or re ceivers in underground mines. 3.5.1 Source and receiver on separate interfaces Let the source S be located on the th interface with 1 < k. The downgoing traveltime is obtained by summing the layer traveltime increments t, i = j,j + 1,. . . , k β 1. 45 Equation (3.7) generalizes to kβi h(O,O)pj, βdown = vi + XP(Thlc,zpk_i,k,z β flk,zPkβ1,k,z) + Yp(flic,zPkβ1,ic,y β βk,Pkβi,k,z) Vklflk,z β β Thj,zPjk,z) + YS(Thj,ZPjk,y β Thj,yPjk,z) (3.14) L vjnj,z Similarly, if the receiver R is located on the 1 interface (1 1 < k), then the upgoing traveltime becomes = β xq(mk,zqk_i,k,z β Thk,zqk_i,k,z) + YQ(nk,zqk__1,k,y β nk,yqk_1,k,z) Vk_lflk,z + XR(flJ,zqjk, β Thj,zqlk,z) + yR(mz,zqzic,y β nz,qzk,z) (3.15) vznl,z The critically refracted traveltime increment is still given by equation (3.9). Summing this expression together with (3.14) and (3.15) yields the total traveltime: 46 h(O,O)pk, h(O,O)q, Tk(XS,YS,ZS,ZR,YR,ZR) = β ___ ___ ___ + x(n1,q11 β nl,,,qlk,z) + yR(nl,zqzk,y β fl,yqlk,z) vlnl,z β XS(flj,zpjk,,, Thj,zPjk,z) + YS(flj,zpjlc,y β 7Zj,yPjk,2) (3 16) vjmi,z where zs = zj(xS,yS) and ZR = zl(rR,yR). This is the proper expression for head wave traveltime when source and receiver are located on different interfaces of the model. It differs significantly from the analogous formula published by Diebold (1987). His expression is actually a special case of the general equation (3.16); in particular, it is only valid if both the source and receiver interfaces are horizontal (:j,z flj,y = Thi,,, = l,y = O flj,z = l,z = 1). The difference between these two formulae is clearly revealed by a detailed examination of a simple two dimensional earth model in Appendix C. The analysis demonstrates that equation (3.16) reduces to the known traveltime solution for this situation, whereas Dieboldβs expression yields an erroneous result. 3.5.2 Arbitrary source and receiver locations A further generalization is obtained by allowing the source and receiver to be located within designated layers. Assume that the source is located in layer j at a vertical depth ds below the the immediately overlying interface (the ih). Similarly, let the receiver be located within layer 1 at a depth dR beneath interface 1. These incremental source and receiver depths must satisfy 0 d h,(as,ys) and 0 dR hl(XR,yR), respectively. The previously developed techniques can be used to derive the head wave traveltime for this situation. Traveltime increments induced by the source layer j on the downward path and the receiver layer 1 on the upward path must be treated separately, because the wave does not propagate across the full thickness of each layer. The result of the analysis is 47 h(O,O)pk, dspk, h(O,O)q, Tk(XS,YS,ZS,XR,YR,ZR) = β ___ __ β + vi vj . vi v + XR(nl,zquc,z β fll,zqzk,z) + yR(nz,zqzk,y flt,yqlk,z) vlnl,z β XS(flj,zpjk,z β Thj,zPjk,z) + yS(flj,zpjk,y β βj,yPjk,z) (317) vjnj,z wherezs = zj(xs,ys)+ds and ZR = zz(XR,yR)+dR. Note that the prior expression (3.16) is recovered in the limit as β 0 and dR β* 0, as expected. As an additional check, examine the case where the source and receiver approach the basal interfaces of their respective layers (interfaces j + 1 and 1 + 1). The layer thicknesses h2(0, 0) and h1(O, 0) at the origin can be rewritten in terms of the thicknesses at the source and receiver positions (h(xs, ys) and hl(XR,yR)) via expression (3.3). Then, by appealing to Snellβs law at each basal interface, equation (3.17) is recast as hi(0, O)Pilc,z [ha(xs, ys) β dsjpk, Tk(XS,YS,ZS,XR,YR,ZR) = . + 1=2+1 7)3 β h(O,0)q, β [hz(XR,yR) β vi vj izzj+1 + xRzl,zqzl,k,, β flli,zqZ-f-1,k,z) + YR(fl+1,zql+1,k,y β vl+lnz+1,z β X8(flj+1,zpj+1,k,z β j+1,zPj+1,k,z) + YS(flj+I,zpj+1,k,y β Thj+1,yPj+1,k,z) . (3.18) vj+lnj+1,z 48 Comparing equations (3.18) and (3.16) indicates that the proper form of the traveltirne expression is obtained in the limit as ds β ha(xs, Ys) and dR β* hl(XR, YR). An additional benefit accrues from separating the downward and upward sums in the traveltime formulae: asymmetric wave propagation paths can be treated. An asymmetric raypath is defined as one where the mode of dowr&going wave propagation in the ith layer differs from the mode of upgoing wave propagation across the same layer. Strictly, differ ent symbols should be used to designate the wave speeds within layer i in the downward and upward sums (e.g., and vβ for the velocities of the dowugoing and upgoing waves, respectively). However, this complication is avoided for the time being in order to main tain notational simplicity. The velocity v appearing in each sum is simply interpreted as the propagation speed of the the desired mode (P or S) across layer i. Equations (3.17) and (3.18) then constitute general formulae for point-to-point traveltimes of head waves propagating in a three dimensional layered earth model. 3.5.3 Arbitrary reference points for layer thickness Individual layer thicknesses enter the traveltime expressions evaluated at the coordinate origin 0. An alternate form of the traveltime equation is characterized by layer thicknesses specified below the source and the receiver. This variant is particularly suitable for the time term, delay time, and reciprocal time inversion methods. Hence, the previous derivation is now modified to incorporate layer thicknesses prescribed at arbitrary reference locations; these points can then be specialized to the source and receiver positions. The resulting traveltime expression forms the point of departure for a three dimensional extension of the aforementioned inversion techniques. The depth of the i interface, referred to an arbitrary location A with horizontal coor dinates (XA, VA), is given by ri,y1 zj(x,y) = zj(XA,yA) β (x β xA)[β--βj β (ii β YA)[zj. (3.19) n:,z ni,z 49 Consider the downgoing raypath first. The previous expression for the traveltime incre ment Lj induced by wave propagation across the th layer is trivially modified to Pilc,z r 1 Pik,y 1 1 Pik,zt = ___ _ β ZA) β (x XA)j + β YA) β (Yi β YA)j + ββ(z+ β z). vi vi vi The vertical (z) coordinates of the ray intersection points are now expressed in terms of the corresponding horizontal coordinates by using equation (3.19): IThizi Ifliy zj β z = h(xA,yA) + (x β X4)Lj + (y β lil,z n:,z lmi+1,zl ___ ___ β (x+ β XA) β (y+i β YA) I Iβ Lflj+1,zJ Lfli+1,zJ where h(xA,yA) zj+1(XA,yA) β zzfrA,YA) is the vertical thickness of layer i below the reference point A. Substituting this expression into the equation for L yields = + (x+i β XA) βi+1,zPik,z β fli+1,zPik,z + (Yi-i-1 β YA) i+1,zPik,y β i+1,yPk,z Vi β (x β XA) i,zPik,z β i,zPik,z β (Vi β VA) i,zPiIc,y i,yPik,z ni,z vi ni,z vi The traveltime increments t are now summed over all layers in the overburden. Applying Snellβs law at each plane interface results in kβi sβ., ___ ___ ___ __ 1down = Vi i:=i 50 + (3p β XA)(Thk,zPk_1,k,, β flk,zPkβi,k,z) + (Yp YA)(βk,zPkβ1,k,y β Thk,yPkβi,k,z) Vk_lflk,z β (XS β XA)(fll,zPlk,z β Th1,xPlk,z) β (Ys YA)(Th1,zPlk,c β ThI,zPlk,z) (3.20) Vi fli ,z A similiar analysis yields the upgoing traveltime. Layer thicknesses are now referred to a different arbitrary position B with horizontal coordinates (XB, YB): kβi β ___ ___ ___ __ LUj,β vi β (aQ β XB)(fl,zqki,k,z Thk,zqlc_1,k,z) β (yΓ§ β YB)(flk,zqk_i,k,y β nk,qk_1,k,z) Vk..4flk,z + (XR β XB)(nl,zqik,z β nl,zqlk,z) + (YR β YB)fr1i,Zq1k,, β fll,yqik,z) β’ (3.21) Vifll,z Finally, the critically refracted traveltime increment is given by a simple alteration to the prior equation (3.9): [(XQ β XB) β (xp β XA)i(fl.k,zPkk β Thk,Pkk,z) = Vkflk,z + β YB) β (Yp β YA)jfrβic,zPkk,y β flk,yPklc,z) Vkflk,z β (XA β XB)(Thk,zPkk,z β flk,zPkk,z) + (YA β YB)frk,zPkk, flk,yPkk,z) (3.22) Vkflk,z 51 The total head wave traveltime is now obtained in the usual manner by summing the contributions along all three raypath portions. Adding equations (3.20), (3.21), and (3.22) gives kβi h(zA,yA)pk, β h:(XB,YB)q:k,z Tk(XS,YS,XR,YR) = 2=1 + (XR β XB)(ni,zqik,z β ni,zqlk,z) + (YR β YB)(Thi,zqik,y β fli,yqik,z) vini,z (z5 β XA)(fli,zPlk,z β fli,zPik,z) + (Us β YA)(fli,zpik,i, fli,yPlk,z) vlnl,z + (ZB β 2A)(flk,zPkk,z β flk,zPklc,z) + (YB β YA)(flk,zPkk,i, β flk,yPkk,z) Vkflk,z + F(xp β XA, YP β YA, XQ β B, YQ β YB), where the function F has been previously defined. Once again, applying Snellβs law at the critically refracting horizon demonstrates that F vanishes identically. Furthermore, it is common (but not mandatory) practice to select the reference points A and B to be coincident with the source S and receiver R, respectively. Thus (ZA, YA) = (xs, Ys) and (XB, YB) = (XR, YR). With these substitutions, the above expression simplifies dramatically to kβi h(xs,ys)pk, β Tk(XS,YS,XR,YR) = + (XR β XS)(Thk,zPkk,z β Thk,zpkk,z) + (YR β YS)(flk,zPkk,y β flk,ypkk,z) . (3.23) Vkflk,2 The two dimensional form of equation (3.23) has some similarity to an analogous expres sion published by Ocola (1972), in that individual layer thicknesses are measured vertically 52 below the source and receiver. However, it is quite different from Palmerβs (1980) 2D tray eltime equation, which currently forms the theoretical basis of the generalized reciprocal method (GRM) of seismic refraction interpretation. Following Ewing et al. (1939), Palmer measures thickness perpendicular to the base of each layer, and thus arrives at a different expression. Nevertheless, the GRM can be formulated on the basis of the two dimensional version of equation (3.23); a distinct advantage is then gained in constructing the refractor depth proffle (see Chapters 5 and 6). - 3.5.4 Reflection traveltime The previous analysis has concentrated on. waves that are critically refracted at the kth subsurface horizon. However, it is straightforward to demonstrate that the traveltime formulae also apply to waves that are reflected from interface h of the model. Snellβs law of reflection for this situation is illustrated in Figure 3.4 and is compactly expressed as k X Pkβ1,k β k X qkβI,k d β u Vkl Note that equation (3.24) allows for a possible mode conversion upon reflection, i.e. v1 is not necessarily equal to v_1. For a wave reflected at interface k, points P and Q are coincident; there is no intervening critically refracted raypath segment. The total traveltime is obtained by simply adding the downgoing and components. Hence, setting (XQ, yq) = (xp, yp) and summing equations (3.14) and (3.15) yields 53 Fig. 3.4. Reflection of a raypath at interface k of the earth model. If a mode conversion occurs, then downgoing velocity vLi differs from upgoing velocity v_1. kβi kβi \β h(O,O)pj \\ ___ ___ ___ __ 1dow+1up = β _ t1 + 2.R(fl,zqZk, β nz,zqk,z) + yR(1,zqzk,y β nz,yqk,z) XS(flj,zpjk,z β 17βj,zPjk,z) + YS(flj,zpjk,y β Thj,yPjk,z) + G(xp,yp), (3.25) where the quantity Gfrp, yp) depends on the horizontal coordinates of the reflection point Pk-1,k qk-1,k k and is given by 54 G(xp,yp) = xp 1 1 d (flk,zpkβi,k,z β β1c,zPkβi,k,z) (nk,zqk_i,k,z β nk,zqk_1,k,z) k,z Vk_i kβi yp 1 1 + d (flk,zpkβ1,k,y β k,yPkβi,k,z) β u (nk,zqk_i,k,y β nk,yqk_i,k,z) k,z 13kl Distinct downgoing and upgoing layer velocities are explicitly incorporated into the above expressions in order to emphasize the possibility of asymmetric mode-converted raypaths. The component form of Snellβs law of reflection implies that G(xp, yp) vanishes. Equation (3.25) is then identical in form to the previous traveltime expression (3.16) that was derived for head waves! 3.6 Rapid traveltime computation In order to compute traveltimes via the above formulae, the unit propagation vectors Pik and qjk overlying the refracting/reflecting interface must be determined. For the reflection problem, this set of vectors depends on both the offset distance X and the azimuth angle 1J1 of the receiver relative to the source. However, the critical refraction problem is qualitatively different; the propagation vectors depend only on the azimuth βIβ. This particular feature can be exploited to yield a rapid computational procedure for head wave traveltimes. Since the propagation vectors depend on the recording azimuth, they should be written as Pk(βP) and qik(), although the explicit dependence on βP is often suppressed for notational convenience. The functional form of this dependence is not known. However, with a minimal amount of raytracing, it is possible to numerically generate the function βP(pjk, q) over the full range of possible recording azimuths (2z radians). Inversion of this function then yields the propagation vectors for a prescribed value of the source-receiver azimuth angle. This technique is discussed in both general and mathematical terms in the following two sections. 55 3.6.1 General description The computational algorithm is based on the close relationship between a critically reflected and critically refracted raypath. For a given source-receiver azimuth angle, both raypath types possess the same unit propagation vectors pk(P) and q,j(J!) for i < k. Consider the following six-step calculation procedure, with reference to the critically re flected/refracted raypath segments depicted in Figure 3.5: (i) Select a point P on the critically refracting interface. Since the location of P is arbitrary, it is convenient to position it beneath the coordinate origin 0. (ii) Choose a critically refracted propagation direction Pkk through point P. The orientation of this vector within the plane of interface k is defined by an angle x measured from an arbitrary reference line. (iii) Forward raytrace from P along the set of upward unit propagation vectors (i = k β 1, k β 2,.. . , 2, 1) to establish the position of a receiving point R on the surface. The departing vector qkl,k is oriented at the critical angle k sinβ(vkl/vk) relative to the interface normal and is contained in the plane defined by Pkk and The appropriate three dimensional form of Snellβs law is applied at all interfaces intervening between the refractor and the surface. (iv) Reverse raytrace from P along the set of downward unit propagation vectors Pik (i = k β 1, k β 2, . . . ,2,1) to establish the position of a source point S on the surface. The incidence angle of the arriving vector pk1,k (in the plane defined by Pkk and nk) also equals the critical angle. A three dimensional βbackward propagatingβ form of Snellβs law is applied at all refracting interfaces. (v) Calculate the azimuth angle of R with respect to S. (vi) Increment angle x by a small amount and repeat steps (ii) through (v). Stop after x has been incremented by a total of ir radians. 56 . . . . Fig. 3.5. Schematic representation of a raypath critically reflected/refracted at point P on interface k. S and R denote a surface source and receiver, respectively. Pk--i,k qk-1,k 57 This procedure numerically defines a function = f(x) over an interval {xo, xo + it]. It is not necessary to perform raytracing to determine the azimuth angle for x E [xo + it, xo + 2ir]. Rather, values are easily generated from the symmetry relation βI = f( β it) + it. This symmetry condition arises from raypath reciprocity: reversing the direction of the critically refracted propagation vector Pkk merely interchanges the positions of the source and receiver on the surface. If the propagation vector Pkk makes one complete rotation on the critically refracting interface (i.e., x increments by 2ir radians), points R and S make one complete closed circuit on the surface (i.e.,, βI also increments by 2w). This functional dependence is designated = (x). The inverse function x = gβ(βP) can then be used to determine the appropriate value of x for a specified source-receiver azimuth angle. Finally, unit propagation vectors Pik and qik corresponding to this value of x are regenerated via points (iii) and (iv) above. Note that the function = g(x) needs to be calculated only once. All recording azimuths contained in the data acquisition geometry are treated by this same function. The numerical method for inverting g(x) is briefly stated in the next section. 3.6.2 Raytracing technique The pertinent mathematical details of the aforementioned computation procedure are now described. A distinction must be made between the downgoing and upgoing portions of the critically reflected raypath. Thus, superscripts βdβ and βuβ refer to these raypath segments, respectively. The Cartesian components of the critically refracted propagation vector Pkk are chosen to be Pkk,z = cosXcoscbkcosOk β sinxsin, (3.26a) = cos x cos sin 0k + sin x cos 6k, (3.26b) Pkk,z = β cossinq. (3.26c) 58 It is straightforward to verify that Pkk 1 and Pkk = 0. Hence, Pkk is a unit vector perpendicular to the interface normal kβ’ Its orientation within the plane of interface k is determined by the angle x If x 0, then Pkk points in the maximum updip direction along the critically refracting horizon. Increasing the angle x rotates Pkk in a clockwise sense when viewed from the coordinate origin 0. Consider the downgoing portion of the raypath first (segment SF). Shah (1973) uses the following three dimensional form of Snellβs law for the ray refracted at interface i: d V1 d Pa βPβi,k + (cos ii β β cos i )n1. Vj_1 V_ Hence, the transmitted propagation vector is a linear combination of the incident propaga tion vector and the interface normal vector. Solving for Piβi,k yields Piβ1,k = + (cos ,.d β cos vd)n1. (3.27) Vt Vt This form gives the incident vector in terms of the transmitted vector and the interface normal. The cosines of the incidence angles can be obtained by the following four-point prescription: (i) cos = Pa . n, (ii) sin d vβl β cos2 vd, (iii) sin = (v_1/ )sin .A, (iv) cos /1 β sin2 Similarly, the upgoing ray (segment FR) at interface i satisfies Snellβs law in the form ii U U qiβ1,k = β--qik + (cosv ββββcos/1 )n, (3.28) Vt Vt where the cosines of the incidence angles are evaluated via: (i) cos cia n, (ii) sin i .,/Tβ cos2 (iii) sin U = (v_1/ ) sin jtL, (iv) cos = β\/1 β 2 yUβ’ Note that a negative sign is used to calculate cos v since the angle v exceeds 7r/2 radians. Starting with pj (= kk) given by expression (3.26), equations (3.27) and (3.28) are evaluated recursively for i = k, k β 1, k β 2,. . . , 3, 2 to generate the unit propagation direction vectors in all layers overlying the critically refracting interface. 59 Intersection points of the critically reflected raypath with the plane interfaces of the model are easily determined after the propagation directions are known. For the downward path segment within the iβ layer, these are related via d dr2 = β 3.29 Hence, r can be calculated if r+1 is already known and the path length l can be found. Taking the dot product with n: and using r n1 = dj (see equation (3.1)) yields 1d β (r+i . n) β d β (Pikfli) Similarly, on the upward propagating raypath segment within the same layer: r = r1 + (3.30) with β (r1.n) (q . n) Thus, starting with r = OI+Oj+zk(O, 0)k, equations (3.29) and (3.30) can be evaluated recursively with i = k β 1, k β 2,. . . , 2, 1 to yield the ray intersection points r and r on the surface. Finally, the azimuth angle Eβ of r (the receiver point) with respect to r (the source point) is calculated by = tan1 β . (3.31) β iJ is obviously a function of the reference angle x: βIβ = g(). In principle, the inverse function x = g1() can be determined. In practice, numerical calculation of the com plete inverse function g1 is unecessary. Rather, inverse interpolation between neighboring 60 computed values of βIβ is used to locate the particular value corresponding to a specified source-receiver azimuth angle βP. Linear interpolation is adequate, provided that the sam pling in x is sufficiently fine. The value is then used to regenerate the complete set of propagation vectors Pik and q needed to evaluate the head wave traveltirne formula. A comparison with an example presented by Chander (1977b) validates the accuracy of the calculation. He considers a four layer model defined by the parameters ci = 0Β°, = 00, zi(0,0) = 0 m, = 4000 m/s, = 0.81Β°, 02 = 225Β°, z2(0,0) = 50 m, v2 = 4500 m/s, q53 = 4.04Β°, 03 = 225Β°, Z3(0,0) = 150 m, 1)3 = 5000 m/s, = 11.98Β°, 04 = 225Β°, z4(0,0) = 350 m, βV4 = 5600 rn/s. Using the above procedure, unit propagation vectors corresponding to the raypath critically refracted on interface 4 are determined for a recording azimuth βP = 33.69Β° (the sampling interval used for x is 0.25Β°). In particular, the downgoing propagation vector at the source βS P14(33.69) = (0.523126)i + (0.322360)j + (0.788938)k. Chander (1977b) iteratively solves two coupled nonlinear equations for the raypath takeoff angles at the source, given the above azimuth angle to the receiver. The orientation of the takeoff vector P14 can be described by two angles in a manner analogous to the interface normals (see Figure 3.1): P14 = (sin cos 714)i + (sin cr4 SLfl y,4)j + (cosa,4)k. (0Β° <a14 < 180Β°) is a polar angle measured from the +z axis and (0Β° -y, <360Β°) is an azimuthal angle measured from the +x axis. Then, the above results yield = 37.91Β° and 31.64Β°, in agreement with Chanderβs (1977b) solution for the same two angles. 61 3.7 Forward modeling examples Two examples presented in this section illustrate the utility, as well as some of the limitations, of the head wave traveltime formulae for forward modeling applications. The recording geometry in the first example is the common reversed proffle; all sources and receivers are located on the surface. The second example considers a typical offset VSP geometry (surface source and downhole receivers). 3.7.1 Profile geometry Equation (3.12) expresses the head wave traveltime in terms of the horizontal offset distance between a surface source and a surface receiver. It is written in condensed form as Tk(xs,ys,X,J!) β mjJ!)X + bk(as,ys,1), (3.32) where the definitions of the slope mk(βI!) and intercept bkQes, ys, are obvious. This ex pression is a three dimensional extension of the slope/intercept formulae that are commonly used to describe head wave traveltimes. Equation (3.32) is evaluated for a shallow three-layer model defined by the parameters qSi 00, 6 = 00, zi(0,O) = 0 m, = 1000 m/s, = 8Β°, 62 = 0Β°, z2(0,0) = 5 m, V2 = 1800 m/s, 63ZZ600, z3(0,0)=12m, v3=3200m/s, Figure 3.6 displays the head wave arrival times observed by a set of four reversed refraction profiles. Direct wave traveltimes are also included in each plot. Shots are located at each end of the recording spreads and the maximum source-receiver offset distance is 50 m. The profile pairs are all centered on the coordinate origin and are oriented in the N-S, NE-SW, E-W, and SE-NW directions (panels a through d, respectively). (m s) (m s) (m s) (m s) N) - N C), 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Co q q C C ) β t . CD N ) 0 p, C D C ; (tp ,β’ β’ β . β b D c- f- p [ 0 a C D L II C β β C n β β s CD : .: CO β C D C β β 4 CD C D C I) C H 63 The straight line arrival time curves plotted in each panel convey an impression of a two dimensional earth model. The full three dlimensionality of the subsurface is only appreciated by comparing traveltime curves recorded along several separate azimuths. Note that the reciprocal times (the shot-to-shot traveltimes) for forward and reverse arrivals in all panels agree, as expected. Refraction traveltime curves are extended to zero offset distance, even though head waves do not exist in the precritical offset zone. Since the traveltime computation method does not locate the critical offset distance, equation (3.32) is evaluated over the full offset range covered by the receiver array. If traveltime analysis is concerned solely with first arrivals, then these calculated nonphysical traveltimes do not pose any problems, because they are always associated with later arrivals. However, precritical offset arrivals do have interpretive significance (Ackermann et al., 1986) and thus their inclusion in the current algorithm may be useful for some studies. Note that the intercept time in equation (3.32) depends on the recording azimuth angle 1β (through the propagation vectors Pik and q) in addition to the source coordinates (xs, ys). This unusual feature appears to be peculiar to three dimensional, multilayered earth models. Intercept time is obviously independent of proffle azimuth for all one dimensional models. For two dimensional earth models (recorded normal to strike), the identity of intercept times observed on split spread proliles is an interpretive rule (Johnson, 1976; Merrick et al., 1978; Ackermann et al., 1986; Briickl, 1987). Finally, in the case of the simplest three dimensional model consisting of a single layer overlying a halfspace, the intercept time is also independent of recording proffle azimuth (see equation (2.10)). Dependence on the azimuth angle β only arises when multiple layers in three dimensions are analyzed. However, raypath reciprocity requires that the propagation direction vectors satisfy Pk(T! + ir) βqzk(1β), qk(P + ir) = Pik(β1β). (3.33) Substituting these results into the expression for the intercept time yields ys, + 7r) bk(XS, ys, P). (3.34) 64 Hence, the two dimensional interpretive rule also holds for the particular three dimensional model examined here (uniform velocity layers bounded by plane interfaces). However, it is not true that that intercept times recorded on all line profiles emanating from the same shotpoint are identical. For realistic earth models (maximum dip 100) this variation in intercept time with profile azimuth appears to be minute. In the current example, b3(0, 0, βP) varies by only 0.03% as the azimuth β] increases from 00 to 3600. 3.7.2 VSP geometry The final example examines head wave traveltimes recorded in a simulated VSP ex periment. The computation procedure described in section 3.6 is readily generalized to a situation where the source and receiver are located on different interfaces of the model. An azimuth angle function β1 = g() can be calculated for a hypothetical source located on interface j and a hypothetical receiver located on interface 1. In an offset VSP survey, the source is located on the surface (j = 1) and a borehole geophone is lowered continuously down the well. Thus, an azimuth angle function is computed for all interfaces 1, including the critically refracting horizon (1 = k). These functions are then inverted at the known source-receiver azimuth angle, and the resulting propagation vectors are used in equation (3.26) to obtain head wave arrival times at each interface intersection with the well. The traveltime to a geophone located within a layer is then calculated by linearly interpolating times computed at the bounding interfaces. This example considers a four layer earth model defined by the parameters: = 00, = 00, zi(0,0) = 0 m, = 1800 m/s, = 30, = 900, z2(0,0) = 50 m, V2 = 2500 m/s, = 40, = 2700, z3(0,0) = 120 m, v = 3200 rn/s, = 00, 04 = 90Β°, z4(0,0) = 155 m, β04 = 3900 rn/s. 65 240 β 220(I) 200 180 240 ββs 220U) 200 180 240 r- 2201.i) 200 180 0 160 Fig. 3.7. Head wave arrival times recorded in a VSP configuration. Surface source is offset from the well by 500 m to the north, east, and west in panels (a), (b), and (c), respectively. 40 80 120 depth (m) 66 The model is strictly two dimensional; the strike directions of all interfaces are north-south. However, since sources are deployed at various azimuths around the well, the three dimen sional formulae are needed to compute accurate arrival times. Figure 3.7 displays head wave traveltime curves as a function of geophone depth. within a borehole positioned at the co ordinate origin. Surface sources are offset from the well by 500 m to the north, east, and west (panels a, b, and c, respectively). For this 2D model, a source offset 500 m to the south generates traveltimes that are identical to those plotted in panel a. The arrival time of a given head wave decreases as the geophone is lowered in the well, because the receiver approaches the critically refracting interface more closely. As the geophone passes through an interface, the slope of the traveltime curve changes. Expressions for these slopes can be obtained from equations (3.17) or (3.18) above. Each curve terminates at the depth of the critical refractor in the well since the associated head wave is not observable below this level. Note that the critically refracted waves depicted in Figure 3.7 are not necessarily the initial arrivals. Other waves that are neglected in the modeling procedure (e.g., direct waves) may actually arrive first over certain ranges of receiver depth. 3.8 Conclusion The traveltime formulae derived in this chapter are useful for a variety of forward mod eling applications involving head waves propagating in three dimensional, multilayered earth models. Obviously, construction of traveltime curves for a trial earth model can assist in the interpretation of field recorded data. The equations also provide means for evaluating the importance of three dimensional effects on head wave traveltimes in various seismological contexts. For example, Merrick et al. (1978) investigate the hidden layer phenomenon using a 2D earth model. Hunter and Puilan (1990), using a 1D model, compare the sensitivities of vertical and horizontal receiver arrays for discriminating layer velocities. Many investi gators are concerned about the possibility that head waves might constitute first arrivals in crosswell transmission tomography experiments. All of these studies can benefit from a full three dimensional analysis. 67 Several straightforward extensions of the results described in this chapter can enhance the utility of the formulae for forward modeling purposes. Inclusion of multiple reflections within the overburden of the critically refracting interface does not pose any special problems (note that multiple raypath segments along the critical horizon are not allowed; see ervenΒ± and Ravindra, 1971, p. 210). Also, well known tools of asymptotic ray theory can be applied to calculate the particle displacement amplitude and waveform of a head wave in this three dimensional situation. Care must be exercised in treating shear wave propagation, because SV and SR modes are not globally decoupled in 3D. Richards et al. (1991) examine some of these phenomena and propose a particular computing procedure. Interestingly, they claim (but do not prove!) that Dieboldβs (1987) traveltime formula is, in general, correct. Thus, their results should be treated with caution. Finally, since head wave traveltime can be expressed by a simple mathematical for mula, inverse methods designed to recover the earth model parameters from measured data are facilitated. Chapter 4 describes an inversion procedure that exploits the rapid forward modeling capability developed above. Other head wave traveltime inversion procedures are probably possible. In particular, the βslope and interceptβ equation (3.32) may allow the ex tension of the two dimensional methods of Dooley (1952), Adachi (1954), and Johnson (1976) to three dimensional models. Although considerable effort has been devoted to discovering this generalization, success has not yet been achieved. 68 CHAPTER 4 INVERSION OF HEAD WAVE TRAVELTIMES FOR THREE DIMENSIONAL PLANAR STRUCTURE 4.1 Introduction Numerous investigators have studied the inversion of reflection traveltimes for three dimensional subsurface structure (Hubral, 1976; GjΓΈystal and Ursin, 1981; Chiu et aL, 1986; Chiu and Stewart, 1987; Lin, 1989; Phadke and Kanasewich, 1990). However, the analogous situation for refraction traveltime data has not been thoroughly examined. Kanasewich and Chiu (198) present a method for jointly inverting reflection and refraction traveltimes to recover 3D structure. Forward calculation of traveltimes is achieved with the iterative ray bending method of Chander (1977a). This chapter describes an algorithm for inverting head wave arrival times for three di mensional planar structure. The earth model is characterized by a stack of uniform velocity layers bounded by plane, dipping interfaces. The inversion method is iterative; an initial estimate of the model parameters is refined until an acceptable match is obtained between observed and predicted traveltimes. Rapid forward modeling of both traveltimes and tray eltime derivatives is achieved with the numerical technique developed in Chapter 3. A novel feature of the inversion procedure is the inclusion of constraint information in the form of in equality relations satisfied by the model parameters. Often, a priori geological or geophysical information is available to guide and constrain a traveltime inversion. This is particularly useful for the inversion of head wave arrival times because the problem can be very ifi-posed and admit numerous solutions. After a discussion of the mathematical basis of the inversion technique, the algorithm is tested on both synthetic and field acquired traveltime data. Inversion of refraction data acquired in the Peace River Arch region of northern Alberta indicates that the algorithm can be a valuable tool for the analysis of traveltimes recorded in a broadside configuration. 69 4.2 Inversion mathematics 4.2.1 General theory Let the observed arrival times for the experiment be organized into an M-dimensional column vector t0. The earth model is characterized by a finite set of scalar parameters m, j = 1,2,. .. , N. These are organized into an N-dimensional column vector m. Predicted head wave traveltimes generated by this model are designated by the M-dimensional vector td(m). Then, a first order Taylor series expansion of the observed data about the particular model m yields the expression A(mβ) zmββ = t(m), (4.1) where zt(mβ) t0,3 β t,,.d(m) is the data discrepancy vector, m1 mt1 β mβ is the parameter update vector, and the elements of the M x N sensitivity matrix A(m?z) are given by [A(mβ)J = 2 (4.2) 23 8m In these and subsequent expressions, the superscript n denotes the iteration index. In many crustal seismic reflection and refraction experiments, the system (4.1) is overde termined, underconstrained, and inconsistent. A popular solution technique for mβ1 is the damped least squares method (Braile, 1973; Kanasewich and Chiu, 1985; Chiu et al., 1986; Chin and Stewart, 1987; Phadke and Kanasewich, 1990). The updated parameter vector is then obtained via mβ1 = mβ + Zmfl+l. Iterations continue until an acceptable fit to the observed traveltime data is achieved. This strategy is termed βcreepingβ (Scales et at, 1990) because the final solution is obtained by the addition of (possibly many) small perturbations to an initial guess. An alternative approach is used in this study to obtain the improved model parameter vector. Substituting the definition of the parameter update vector into (4.1) and rearranging 70 terms gives At(m) + A(mTh)m. (4.3) The right hand side of this expression consists of known quantities. Hence, equation (4.3) can be solved directly for the new model parameter vector mβ1. Scales et al. (1990) refer to this technique as βjumpingβ because, in the absence of additional constraints, the size of the model change between iterations m β JJ is not restricted to be small. Since the inversion is formulated in terms of the model itself, rather than a model perturbation, the jumping strategy facilitates the incorporation of constraint information into the algorithm. Constraints are mathematically expressed in the form of inequality relations satisfied by the model parameter vector: m_ < m1 m, (4.4a) where the vectors m_ and m+ are lower and upper bounds on the model, respectively. These bounds arise from a priori geological or geophysical knowledge (or assumptions) about the earth model. For example, the model parameters in this study are strictly nonegative. Thus, if the lower bound is set equal to 0, negative values are excluded on each iteration of the inversion procedure. This is required for meaningful forward modeling of traveltimes. The inequality bounds in (4.4a) can also be used to severely restrict (or even eliminate) the variation of a certain parameter on succesive iterations of the inversion. In this situation, upper and lower bounds are narrowly established about an accurately known (or preferred) value for the particular parameter. The inversion algorithm is also stabilized by limiting the size of the model change between iterations. Hence, if Sm is a vector of upper bounds on the parameter increments, the updated parameter vector must satisfy the additional inequality constraints m β Sm < m < m + Sm. (4.4b) These constraints fulfill the same regularizing role as the damping parameter in a damped least squares solution of the original equation (4.1). If a reasonably good initial estimate 71 m0 for the model is available, and if II Sm II is sufficiently small, then the constraints (4.4b) assure that the algorithm iterates toward a solution in the neighborhood of m0, rather than jumping to a remote region of model space. The observed traveltimes on the right hand side of the linearized data equation (4.3) are contaminated with random picking errors. Hence, an exact solution may not exist. A robust solution to this inconsistent system can be obtained by minimizing the l norm of the misfit. Linear programming provides a convenient solution method because the model parameters are intrinsically non-negative and are constrained by inequalities (4.4a,b). In order to pose the problem in the context of linear programming, an M-dimensional residual vector r is introduced into equation (4.3) as follows: A(m) mβ1 + r = zt(m) + A(mTh)m. (4.4c) The elements of r constitute additional unknown variables that must be solved for. The problem now consists of determining the model parameter vector m1 and the residual vector r that simultaneously satisfy the inequality constraints (4.4a,b), the equality con straints (4.4c), and that minimize the 1 norm of the residual r = (4.4d) A standard linear programming routine is used here to solve the constrained optimization problem specified by equations (4.4a,b,c,d). After an improved model parameter vector mβ is obtained, the lj norm of the misfit between observed and predicted traveltimes tob3 β t,.d(mββ) is computed. Iterations cease when the misfit reaches some acceptable level, or exhibits negligable change on successive iterations. The elements of the residual vector r play no further role after equations (4.4) are solved, and are discarded. However, the 11 norm of the residual (4.4d) is monitored on each iteration in order to assess how closely the linearized data equations have been fit. 72 Finally, the flexibility of the inversion algorithm is enhanced by including variable weight ing of both the traveltime data and the model parameters. Equation (4.3) is modified to [wd A(m) W1] [w m1] = Wd t(mβ) + [wd A(mβ) wiβ] [w mj, where Wd and W are data and parameter weighting matrices, respectively. Currently, these are restricted to be diagonal matrices. Thus, premultiplication of the sensitivity matrix by Wd corresponds to row weighting and postmultiplication by Wβ corresponds to column weighting. The data weighting matrix can be used to emphasize those particular traveltimes judged to be more significant for the inversion. The parameter weighting matrix serves to nondi mensionaiize and normalize the elements of the model parameter vector m1. This is a practical concern in inversion algorithms where the model is characterized by parameters with different physical dimensions and/or widely varying numerical magnitudes. The con ditioning of the sensitivity matrix A(mβ) is improved by column scaling by Wiβ. Thus, numerical roundoff error associated with the linear programming solution is reduced. Suit able units of measure are chosen for the various model parameters; the reciprocals of these scalars form the diagonal elements of Wi,. The inequality constraints (4.4a,b) on the model must also be nondimensionalized in the same manner. If the weighted parameter vector calculated by the inversion algorithm is designated thL+l, then the physical parameters required for forward modeling of traveltimes are obtained via mβ = W;1 thββ1. 4.2.2 Calculation of sensitivities A simple 3D earth model consisting of a single layer overlying a halfspace is characterized by the five-element parameter vector m = [vi,v2,qf,O,hjTwhere v and v2 are P-wave velocities of the two media, and 0 are interface orientation angles, and h is the vertical depth to the refractor below the coordinate origin (Figure 2.1). A closed form expression for head wave traveltime in terms of these parameters is given in Chapter 2 (equation (2.10)). Hence, formulae for the elements of the sensitivity matrix can be derived by straightforward 73 partial differentiation. However, a substantial simplification arises by redefining the set of model parameters as m = [SI, c, b, 8, hJT, where sj = 1/vi is the layer slowness and = sin(vi/v2) is the critical refraction angle. The traveltime equation is then linear in the slowness Si. In terms of these new parameters, head wave traveltime is T = sisin(i β6)X + 2sicosic[hcoscβsinc(xscos6+yssin8)], (4.5) where S sin1 [sin 4. cos(111 β 6)]. z5 and yg are the horizontal coordinates of the source, - and X and βIβ are the offset and azimuth to the receiver, respectively (Figure 2.2a). Differ entiating expression (4.5) yields the traveltime sensitivities = sin(iβ8)X + 2cosic[hcos4._sinc5(xscos8+vssin8)], = sicos(iβS)X β 2sisinic[hcos4.βsin4.(xscos8+ssin8)], 8T βsicos(iβS)cos4.cos(PβO) . r . . i β = X β 2si coszihsmc5 + cos 4.(xs cos9 + y sin 0)i,84. cosS L J ΓT βs cos(i β 6)sin4.sin(4β β6) β = X + 2si cosz sin q5frssm8 β cos 8),88 cosS ΓT - = 2sicoszcosc5. The i1 row of the M x 5 sensitivity matrix A(m) is obtained by evaluating these relations with the current model vector mlZ and with the geometric parameters {x5,yg2,Xj, βIβ} appropriate for the th recorded traveltime. The number of traveltimes M usually exceeds 5, and thus the matrix is overdetermined. However, in some situations, the matrix is also rank deficient. If the refracting horizon is horizontal (4. = 0) then the derivative 8T/88 vanishes, and the fourth column of the sensitivity matrix is identically zero. This will occur, for 74 example, if a one dimensional earth model is used for the initial parameter vector estimate m0. Calculation of head wave traveltime sensitivities for a multilayered earth model is more complicated because a closed form mathematical expression for the arrival time does not exist. However, Chapter 3 describes a rapid computing procedure for obtaining critically refracted traveltimes for arbitrary recording geometries. A combination of analytical and numerical techniques can then be used to calculate the sensitivities. In. this study, sources and receivers are restricted to the surface. Moreover, the orientation angles of the surface (4i and 8j) are assumed to be known. Thus, if there are K interfaces in the earth model, then there are a total of N = 4K β3 model parameters (K layer slownesses, 2K β2 interface orientation angles, and K β 1 layer thicknesses). These parameters are organized into the N-dimensional column vector m = [sl,...,sK, 82,...,K, hl,...,hK_l], where s1 is a layer slowness, qS is an interface dip angle, 8 is an interface azimuth angle, and h1 is a vertical layer thickness measured at the origin. The numbering convention for layers and interfaces is described at the beginning of Chapter 3. A βslope and interceptβ expression for the traveltime of a head wave critically refracted at interface k (2 k K) is given by equation (3.12). In the current notation, it is T = m(βP)X + b(x,y,βJI), (4.6a) where the slope m(P) and intercept time bfrs, ys, IIβ) are m(βTF) = s [cos iJ qlk,z + sin βIβ qik,y β tan q cos(P β 0β) ik,z], (4.6b) 75 ys, = s h (Pik,z β qik,z) β Si [xs (pik,z β qik,z) + ys (p β qik,y) β tan c5i (xs cos O + ys sin 6β) (plk,z β qik,z)]. (4.6c) Note that the slope and intercept simplify if the surface is horizontal (4β, = 0). Equation -(4.6) is not a βclosed formβ expression for head wave traveltime because the unit propagation vectors Pk and characterizing the raypath depend implicitly on both the model m and the recording azimuth βEβ. Nevertheless, if m and are specified, then all of the propagation vectors can be accurately computed. Many of the elements of the sensitivity matrix A(m) can be determined by analysis. The traveltime of a head wave formed on interface k does not depend on the parameters s, 4βz, and 6 for 1> k, or on h1 for 1 k. Hence, these traveltime derivatives are identically zero. Moreover, the sensitivity to layer thickness h1 for 1 < k is derived directly from the above expression for intercept time: ΓT/0h1 = Si (pi β quz). The remaining sensitivities must be evaluated by a finite-difference technique. Let dm1 be a model perturbation vector with zeros in all element positions except the thβ’ Then, the partial derivative of head wave traveltime with respect to parameter m is approximated by the forward finite-difference 07β T(m + dm,) β T(m) (47) β IIdmiII Unit propagation vectors are generated for both the perturbed model m + dm, and the unperturbed model m via the procedure described in Chapter 3. These are then used in formulae (4.6a,b,c) above to calculate traveltimes for each model. Finally, substituting these traveltimes into (4.7) yields the required derivative. The size of the model perturbation dm, is typically about 1% of the value of the associated model parameter. Although a centered finite-difference scheme would yield greater accuracy, the one-sided approximation 76 used here requires less computational effort. Only two traveltime calculations are necessary for each model parameter, instead of three. 4.3 Synthetic examples The two examples discussed in this section are representative of a large number of computational experiments conducted with the inversion algorithm. Both the single-layer and multiple-layer variants of the algorithm are examined. The first example demonstrates that the single-layer version is capable of returning the correct solution under a variety of operating conditions. However, as indicated in the second example, the multilayer version appears to require fairly restrictive constraints in order to iterate toward the correct model. 4.3.1 Single layer The earth model used for the first example consists of a single layer overlying a halfspace and is defined by the parameter values: vj 1500 m/s, = 2500 m/s, = 50, 9 = 450, h = 100 β’ Critically refracted arrival times for two areal recording geometries are generated from equa tion (2.10). Inversions are performed with both accurate and error contaminated traveltimes. Figure 4.la is a plan view of a triangular data acquisition geometry. Thirty receivers are deployed around the perimeter of an equilateral triangle with sides 500 m long. Each side contains 11 receivers separated by 50 m. All receivers record energy from a source that is activated sequentially at the three vertices of the triangle. However, some source-receiver offsets are less than the critical offset distance (equation (2.11)) and thus the receivers do not detect a head wave arrival. These fictitious times are excluded from the inversions. A total of M = 73 uniformly weighted traveltimes are used to recover the N = 5 earth model parameters. 77 I 100 + + 0β + + ++++ + +++ *+ + + + + + + + +* β200 I I I I I β200 0 200 200b ++++ 100-+ 1 + β200 0 200 easting (m) Fig. 4.1. Plan views of two areal recording arrays. Sources are indicated by asterisks and receivers by small crosses. (a) Triangular geometry. (b) Swath geometry. The strike and clip symbol in the upper right corner of each panel refers to a subsurface interface located 100 m below the coordinate origin (large cross). 78 Numerical results from a few typical inversion runs are described here. The iterative inversion procedure is initiated with a one dimensional model given by vi =l000m/s, v2==2000m/s, q=OΒ°, 8=0Β°, h=85m. After eight iterations, an excellent estimate of the true earth model is returned: v = 1500.21 m/s, v = 2499.98 m/s, = 5Β°, 8 = 45Β°, h = 100.02 m. The initial traveltime misfit of 79.6 ms is reduced to 0.1 ms; iterations are terminated when the relative change in the misfit is less than 1%. A wide variety of starting models yields essentially the same final solution, although the number of iterations required for convergence varies. Also, the inversion is stable when the traveltimes are contaminated with small random errors. Random numbers drawn from a uniform probability distribution on Β±4 ms are added to the accurate times, and the algorithm is initiated with the same starting model. After five iterations, the following solution is obtained: vi = 1496.93 m/s, V2 = 2549.46 m/s, = 5.22Β°, 8 = 51.93Β°, h = 99.65 m. Iterations cease when the misfit decreases below 2.3 ms, or one standard deviation of the noise. Note that the overburden velocity vi has been correctly estimated from the refraction data alone! This interesting (and unusual) result is consistent with the theoretical analysis presented in Chapter 2. These results suggest that the triangular recording array is a useful tool for determining three dimensional planar structure. This particular geometry combines an adequate distri bution of offset and azimuth with three reciprocal time pairs. These are favorable attributes for a successful inversion of refraction arrival times via the time term method (Scheidegger and Willmore, 1957; Wilhnore and Bancroft, 1960; Berry and West, 1966; Smith et al., 1966; Reiter, 1970). The triangle was first investigated by Gardner (1939) who demonstrated that 79 it yields an exact solution for the delay times at the three vertices. More recently, the trian gular array has been used for deep crustal seismic exploration in Saskatchewan (Kanasewich and Chiu, 1985) and British Columbia (Zelt et aL, 1990). Other recording geometries are considerably less robust in detecting and resolving the three dimensional dipping structure. Figure 4.lb depicts two parallel line arrays, separated by 200 m, oriented along the strike direction of the subsurface refractor (Nw-SE). Each spread contains 11 geophones (receiver interval 50 m) that record head wave arrivals from a source located at the center of the opposite array. This broadside recording pattern simulates aspects of the βswath geometryβ commonly used for 3D seismic reflection surveys. Using the same starting model, inversion of the 22 broadside traveltimes yields vj = 1285.85 m/s, V2 = 2500.37 m/s, Γ§t = 4.96Β°, 0 = 45.02Β°, h = 80.00 m. The initial traveltime misfit of 66.4 ms is reduced to 0 ms in eight iterations. Although this model generates an exact fit to the data, the overburden velocity v and vertical depth h are incorrect. A similar effect is observed when the inversion is initiated from numerous different starting models. This situation illustrates the classical tradeoff between overburden velocity and refractor depth in seismic refraction interpretation. If additional a priori data are introduced into the inverse problem, then a correct solution is possible. For example, the interface depth may be known from a borehole drilled at the coordinate origin. Constraining the depth to satisfy 99 m < h < 101 m yields the model = 1490.01 m/s, V2 = 2500.31 m/s, Γ§ = 5.12Β°, 6 = 45Β°, h = 99.00 m, which is substantially correct. Alternately, constraining the overburden velocity with the inequalities 1450 rn/s <v1 <1550 rn/s yields the erroneous one dimensional model β 1450 m/s, v2 = 2493.96 m/s, Γ§i5 = 0Β°, 0 0Β°, h = 94.25 m. 80 Evidently, the broadside geometry illustrated in Figure 4.lb allows many solutions to this nonlinear traveltime inverse problem. However, if the inline head wave arrivals recorded by each spread are also included in the inversion (M = 32) then the correct model is recovered in seven iterations: = 1499.85 m/s, V2 = 2499.99 m/s, qS = 5Β°, 6 = 45Β°, h = 99.98 m. Head waves recorded along strike provide excellent control on the velocity of the refracting medium. Equation (2.13) and Figure 2.4 indicate that the measured apparent velocity equals the true velocity in this situation. As expected, the inversion is degraded when the exact traveltimes are contaminated with uniformly distributed random errors on Β±4 ms. The model = 1558.76 m/s, v2 = 2664 m/s, 4 = 4.99Β°, 6 = 42.72Β°, h = 106.86 m, is obtained in seven iterations with a traveltime misfit of 2.4 ms. 4.3.2 Multiple layers The earth model used for the second example consists of two layers overlying a halfspace and is defined by the parameter values v = 1500 m/s, ci = 0Β°, = 0, = 0 m, V2 = 2000 m/s, qf2:= 3Β°, 62 = 180Β°, Z2 = 40 m, V3 = 2500 m/s, = 50, 63 = 45Β°, = 100 m. As indicated previously, the model parameters qi, Oi, and zi are not allowed to vary in the inversion. Hence, there are only N = 9 parameters to estimate. Head wave traveltimes for 81 each critically refracting horizon are generated by evaluating equation (4.6a). The recording geometry used is the triangular array displayed in Figure 4.la. Also, precritical offset arrivals are included in the traveltime dataset (in an actual field experiment, these times could be estimated by extrapolation or phantoming). Hence, there are a total of M = 174 arrival times input to the inversion procedure. If restrictive constraints are imposed on the model parameters, then it is possible to re cover the correct solution with the multilayer inversion algorithm. For example, constraining the velocity vl and the depth Z2 to be equal to the true values allows the algorithm to iterate - to the known solution from a one dimensional starting model. Also, if narrow bounds are placed on the velocities and depths (Β±50 rn/s and Β±5 m about the true values, respectively) then the correct solution is obtained from a nearby initial model. However, a relatively unconstrained inversion invariably yields an erroneous result for this simulated experiment. It is probable that the- objective function (4.4d), which only measures the misfit in the lin earized data equations, has many local minima that preclude convergence to the desired global minimum. 4.4 Field data example 441 Peace River Arch broadside data Deep seismic refraction data were acquired in the Peace River Arch (PRA) region of northern Alberta in 1985. In addition to four inline profiles, two broadside profiles were recorded. Figure 4.2 illustrates the source-receiver geometry for these two proffles. Shot A4 in the west is recorded by the north-south trending line A, and shot B4 in the north is recorded by the east-west trending line B. Offset distances range from 249-313 km on line A, and from 262-344 km on line B. Each receiver array subtends an azimuthal angle of 64Β° relative to its source. The first arrivals at the recording sites are interpreted to be waves that are critically or near-critically refracted at the Moho (Zelt, 1989). Hence, inversion of the first break traveltimes can provide an estimate of the regional depth and dip of the Moho beneath the Peace River Arch. 82 -c C C * Shot 84 300- 200- 100 - 0 β100 β200 Line B Shot A4 Line A β200 β100 0 100 200 easting (km) Fig. 4.2. Broadside recording geometry for the Peace River Arch seismic experiment. Sources are indicated by asterisks and receivers by small crosses. 83 In this example, the broadside arrival times are inverted within the framework of a simple βlayer over a halfspaceβ earth model. No attempt is made to infer a detailed structural picture of the crust and Moho in the Peace River Arch region. Rather, the intent is to recover a large scale 3D model that can be subsequently refined by other traveltime interpretation/inversion methods. Initial estimates of the model parameters are obtained from crustal sections along lines A and B given by Zelt (1989, p. 103 and 112). These sections were derived by interpreting the inline refraction data of the PRA experiment via a trial-and-error forward modeling approach. There are 83 first break picks from line A and 52 first break picks from line B. Thus, M = 135 uniformly weighted traveltimes are input to the iterative inversion procedure. Since the first arrivals recorded at these long offset distances are emergent, the estimated picking error is relatively large (Β±50 ms). 442 Static corrections The effect of variable near surface structure on the head wave arrival times can be reduced be applying static corrections to the first break picks. Static corrections are commonly applied to seismic reflection data for this same purpose. However, there is an important distinction between the corrections for the two types of data. In the refraction case, the static is designed to remove the refraction delay time influence of the near surface structure, and then replace it with the delay time contribution of a constant velocity medium. In the reflection case, the correction pertains to the vertical traveltime through the actual and replacement media. Statics application is an important preprocessing step for the PRA broadside data because the subsequent inversion assumes a very simple earth model. In effect, an attempt is made to βmake the data fit the modelβ more closely. The following development assumes that near surface velocity information is available. For the PRA experiment, this is obtained from well log data along the two proffle lines. An expression for the static correction at each receiver site is derived by approximating the local near surface velocity structure by a one dimensional stack of layers. The th layer is bounded from above and below by plane, horizontal interfaces at depths z, a nd zj, 84 respectively. The P-wave velocity of the layer depends linearly on depth z and is given by v(z) v(O) + a,z. In this situation, the raypath associated with a particular ray parameter p is a circular arc. The upward propagation time across the layer is it(p) = -- cosh (_2_) β cosh1 (β_) , (4.8) a1 pv1 pv where v v(z) = v.(O) + az and v vj(zj+i) = v(O) + a.jzj+1. Symbols v and v denote the velocities at the top and bottom of the i layer, respectively. The horizontal distance accumulated by this wave traversing layer i is IXj(p) = β-- [βi - (pvj)2 - - (PvD2]. (4.9) If the wave is critically refracted on a horizontal interface at greater depth, then the ray parameter p equals 1/va, where v is the critical refraction velocity. The angle i(z) that the raypath makes with the vertical at any depth is then determined by a velocity ratio: 77(z) sin1 [pvi(z)] = sin1[v1(z)/vj. Using this result, equations (4.8) and (4.9) reduce = (1 + cos ) n = [cos β cos (4.lOa, b) a1 (1 + cos j) sin a where ?7 = sin(v/v) and = sin1(v/v). 77 and 77 are the incident angles of the circular raypath segment at the top and bottom of layer i, respectively. The refraction delay time contribution of the th layer is defined as Jt β 85 Substituting from equations (4.10) yields = ---ln , (4.11) a1 where quantities Q and Q are given by β 1 + cos β 1 + cosQiβ . I. sin i exp(cos ) sin exp(cos i) As a check on the correctness of this result, examine the case where the vertical velocity gradient of the layer vanishes (a β* 0). In this situation, the raypath becomes a straight line segment, and the angles m and i both approach the same angle βo = sin [v(o)/v]. Application of LβHopitalβs- rule to equation (4.11) then yields (z+i β z) cos 77o coshmzr= = v(O) v(0) This is the proper expression for the one-way delay time contribution of a layer with thickness h and uniform velocity v(0). The one-way refraction delay time associated with a stack of m horizontal layers is obtained by summing expression (4.11): tT = = (4.12) This formula naturally accommodates any velocity discontinuities at the interfaces (i.e., v_1 v2). However, if the velocity function is continuous between all layers (as in the PRA 86 example), then Q = Q+i and expression (4.12) simplifies to = E--1n . (4.13)a Q+i The static correction applied to the picked arrival times is designed to remove the one- way refraction delay time influence of the layered near surface structure, and add the delay time associated with a uniform replacement medium of equal thickness. If the replacement velocity is designated Vr, then the static shift is = βIr + (ih)cos. (4.14) where = sinβ (v,./v). Although this formula has been derived within the context of upward wave propagation through a set of horizontal layers, it also applies to the downward propagating portion of the total raypath. Hence, it may be used for the computation of either a source site static or a receiver site static. Near-surface velocity information along lines A and B of the PRA experiment are ob tained from models developed by Zelt (1989, p. 40) from well logs. There are eight wells along line A and seven wells along line B. At each weilsite, the velocity model is approx imated by a one dimensional stack of layers with linear velocities, and a static correction is calculated via formula (4.14). The critical and replacement velocities are assumed to be 8.25 km/s and v,. = 6.6 km/s, respectively. These weilsite static corrections are then interpolated/extrapolated to the receiver locations on each line using the following proceβ dure. First, a straight line is fitted to the coordinates of all recording stations and weilsites along a proffle using the York algorithm (York, 1966, 1967, 1969; Williamson, 1968). This is the proper line fitting technique to use in this situation because both the northing and easting coordinates are subject to positioning errors. The York algorithm also automatically returns the coordinates of the perpendicular projections of the points to the fitted straight 87 line. Finally, the weilsite statics are linearly interpolated/extrapolated in these projected co ordinates to the receiver sites. Figure 4.3a displays the computed receiver static corrections for both lines; squares denote the weilsite statics calculated from expression (4.14). The magnitudes of the statics are approximately the same on both lines ( 300 ms), although there is a noticeable increase toward the southern end of line A. The static corrections calculated with the above procedure apply only to receiver sites. However, since shot A4 is located at the end of profile line B, a source static is readily estimated for this shot from the nearest computed receiver correction (β290 ms). A source static for the isolated shot B4 must be assumed; the value adopted here is β300 ms. 4.4.3. Inversion results Initial estimates of the model parameters are inferred from the interpreted crustal sec tions in Zelt (1989, p. 103 and 112). These starting parameter values are v = 6.5 km/s, vZ = 8.25 km/s, g 2Β°, 8 = 270Β°, Ii = 40 km. The same sections are also used to provide lower and upper bounds on the earth model parameters: vj = 6.0 km/s, v = 7.5 km/s, 4 = 0Β°, 8 = 0Β°, h = 38 km, vt = 7.5 km/s, vt = 8.5 km/s, q5 = 15Β°, 6 = 360Β°, h = 42 km. The bounds on the velocities must be transformed to equivalent bounds on the parameters Sj (slowness) and i (critical angle) that are used by the inversion algorithm. After five iterations, the following model is returned: = 6.21 km/s, v2 = 8.50 km/s, = 2.75Β°, 8 = 299.74Β°, h = 42 km. EE 88 o β100 200 300 β0.20 - 0 100 200 300 β0.20 L : 0 100 200 300 0 100 200 300 inflne distance (km) Fig. 4.3. Receiver static functions for lines A and B of the PRA experiment. Squares refer to well locations where near surface velocity logs are available. (a) Statics calculated with v = 8.25 km/s and yr = 6.6 km/s. (b) Statics calculated with v = 8.5 km/s and v1. = 6.2 km/s. 89 The initial traveltime misfit of 658 ms is reduced to 172 ms; iterations cease when the relative change in the misfit is less than 1%. The velocities vi and V2 of this solution differ from the replacement and critical velocities previously used for calculating the static corrections. Hence, the statics are recomputed with V7 = 6.2 km/s and v, = 8.5 km/s and applied again to the picked arrival times. Figure 4.3b ifiustrates this second set of receiver static corrections. New source statics are assumed to be β210 ms for shot A4 and β250 ms for shot B4. Then, initiating the inversion algorithm with the same starting model yields = 6.17 km/s, V2 = 8.50 km/s, S = 300.86Β°, h = 42 km, in five iterations with a misfit of 171 ms. Although this model is not significantly different from the previous one, it is consistent with the assumptions used for calculating statics. A comparison between the observed arrival times (after static corrections) and the travel- times predicted by the model produced by the inversion is displayed in Figure 4.4. Evidently, the simple five parameter earth model provides an adequate explanation for the gross char acter of the broadside arrival time curves. Smail scale variations in the predicted times (solid curves) are due strictly to recording geometry irregularities, rather than any subsur face structural complications. However, a large component of the total traveltime misfit must be attributed to structure or velocity variations that are not modeled in the inversion procedure. This misfit is too large to be accounted for by random picking errors alone. For example, the predicted times are systematically greater than the observed times throughout the central portion of line B. This suggests that the Moho north of line B is not adequately represented by a plane interface bounded by uniform velocities. Vertical depths to the Moho calculated from the inversion results are posted on a plan view of the PRA recording geometry in Figure 4.5. A depth trend is readily apparent, although the depth at the southern end of line A is probably too large. A southeastward dip of the Moho is suggested by interpretations of the inline refraction data on lines A and B, but not on the other lines of the PRA experiment (Zelt, 1989). Also, the velocities recovered by the 3D inversion are broadly consistent with those obtained by Zelt (1989). Of course, inflne distance (km) Fig. 4.4. Comparison between predicted traveltimes (solid lines) and observed traveltimes (triangles) on broadside lines A and B of the PRA experiment. 90 Line A0 0 48 46 44 42 40 38 50 48 46 44 42 40 38 I I I I I 0 100 200 300 Line B I I I I I I I 0 100 200 300 91 Fig. 4.5. Vertical depths (km) to the Moho in the PRA region inferred fom the inversion results. 300 200- 100 - 0- β100 - β200 0β -c -I-? 0 C 42. * :34.6 42.0 47.2 37.7 47,8 I I I β200 β100 0 100 200 easting (km) 92 precise agreement cannot be expected. The mean crustal velocity is lower (6.17 km/s vs. 6.6 km/s) and the sub-Moho velocity is higher (8.5 km/s vs. β-.. 8.2 km/s). However, Zelt inferred a P velocity of 8.4 km/s along the northern half of line A and the eastern quarter of line B. An independent check on the validity of the 3D inversion is provided by a method proposed by Zelt (1989). Using a simple two dimensional earth model, each arrival time recorded along a broadside proffle can be inverted for an estimate of Moho depth beneath the associated receiver site. A depth proffle for the Moho can then be constructed by plotting these depth estimates side-by-side. Zelt (1989) accomplished this inversion using ray tracing techniques. However, as the following derivation demonstrates, a straightforward mathematical solution to the problem is possible. In the two dimensional situation, head wave traveltime is T(X) sin(i+ + 2h(x) cos βp cos (4.15) where i = sin1(vl/v2) is the critical angle and h(xs) is the vertical thickness of the layer at the source location. The interface dip angle βp may be positive, zero, or negative. Substituting h(x) = h(O) + xstan βp into this expression yields T(X) = sin(i+βp) + 2h(O) C:: βp cos c + 2 sinβp cos h(O) is the layer thickness at the coordinate origin. This relation is rewritten as a quadratic form in cos Acos2βp + Bcosβp + C = 0. (4.16) 93 The three coefficients depend on earth model parameters {vl,v2, h(O)} and measured data {xs,X,T}: A = 1 + 4h(o)cosi [ni + [1 + (i + )jJ v1T F 2h(O) 1B = β2 βiβ Lsmn z + x cos 2 viTi 2x 2 c= β 1+_ cosic. If values for v, V2, and h(O) are known (or assumed), then these coefficients can be evaluated numerically. Solution of equation (4.16) for cos is via the quadratic formula: -B+/B24AC cos(p = 2A β (4.17) where the positive root is chosen by analyzing the form of the right hand side as X β* +00. There is still a two-fold ambiguity in determining the dip angle Γ§o from cos Γ§o. This ambiguity is easily resolved by ensuring that the traveltime predicted by formula (4.15) agrees with the known traveltime T. Finally, once the dip angle is determined, the vertical depth to the refractor at any inline position z can be calculated via h(z) = h(O) + tan p. Figure 4.6 compares Moho depths calculated via the above method with those inferred from the parameters recovered by the 3D inversion. Vertical depths beneath the recording stations on each line are plotted. The coefficients in the quadratic (4.16) are evaluated with the parameters obtained from the 3D inversion procedure (vi = 6.17 km/s, V2 = 8.50 km/s, h(0) = 37.6 km for shot A4, and h(0) = 34.8 km for shot B4). Short wavelength variations in the computed curves are artifacts of the traveltime picking errors, and should be ignored. There is close agreement in the Moho depth trends calculated by these two completely 94 different techniques, especially along line A. The larger departure of the two depth curves along line B suggests structural complexity north of this line in the Peace River Arch region. Although the two dimensional inversion method is simple and appears to yield reason able depth estimates, it requires the assumption of numerical values for the three unknown model parameters {v, V2, h(O)}. In contrast, the 3D inversion technique yields simultaneous estimates of all of the relevant earth model parameters. Moreover, since it accomplishes a joint inversion of all of the error contaminated traveltime data, it is more robust than the 2D method. 4L1neA O 100 200 300 0 100 200 300 inβine distance (km) Fig. 4.6. Comparison of vertical depths to the Moho beneath the receivers of lines A and B calculated by two different methods. Smooth curves in each panel are depths inferred from 3D inversion results. Jagged curves are obtained from a simple 2D inversion method. 4.5 Conclusion The inversion of head wave arrival times for three dimensional planar structure is posed as a constrained, parameter optimization problem. The iterative inversion algorithm de 95 scribed here has the abffity to converge to a realistic solution provided that (i) the data ac quisition geometry is adequate, and (ii) sufficient a priori constraints are available. However, precise definitions of βadequateβ and βsufficientβ in this context are not known. Nevertheless, the inversion procedure provides a useful tool for examining these phenomena. In particular, the single-layer variant of the algorithm can be used as an aid in designing recording geome tries for detecting and resolving three dimensional dipping structure. Also, investigation of both synthetic and field recorded datasets indicates that it can be successfully applied to the inversion of broadside refraction data. Currently, there is a lack of techniques for effective interpretation of such data. Although the multiple-layer version of the algorithm exhibits a greater tendency to converge to an erroneous result, several successful inversions have been achieved with the inclusion of sufficient constraints. Generalization of these algorithms to three dimensional recording geometries would be a useful extension. As indicated in Chapter 3, the forward modeling procedure can easily accommodate buried sources and/or receivers. Moreover, well log information regarding interface depths and layer velocities would provide the necessary constraint information for the multilayer algorithm. 96 CHAPTER 5 POINT DEPTH AND DIP ESTIMATES FROM REFRACTION TRAVELTIMES 5.1 Introduction The Generalized Reciprocal Method (GRM) is a technique for delineating an undulating subsurface interface from refraction traveltime data recorded by inline forward and reverse proffles. After the refracted arrivals from a particular horizon are identified, simple linear combinations of the picked times are used to form two analysis functions: (i) the velocity analysis function, and (ii) the generalized time-depth function. Each is plotted with respect to receiver position along the profile line. The slope of the velocity analysis function is a measure of the P-wave velocity of the refracting layer. The generalized time-depth function is a measure of the depth to the critically refracting interface, expressed in units of time. Hence, it is analogous to a reflection event on a seismic reflection time section. Conversion of the time-depth function to actual depth values requires velocity information for both the overburden and the refracting layer. The GRM was developed by Palmer (1980, 1981) as an extension of the conventional reciprocal time method for interpreting refraction arrival times (Hawkins, 1961). Although it has been successfully applied to the problem of mapping undulating refractors, the math ematical formulation of the method is based on a two dimensional earth model with plane interfaces. Palmer adopts a model parameterization originally proposed by Ewing et al. (1939): the thickness of a layer is measured normal to its basal interface. Moreover, Palmer (1980, p. 3) claims that an alternative parameterization characterized by vertical layer thick nesses is βnot suitableβ for the derivation of the ORM analysis functions. This statement is incorrect. All of the GRM analysis tools can be rigorously derived from the two dimensional version of equation (3.23) for head wave traveltime. This development is the subject of the first portion of this chapter. The use of vertical layer thicknesses results in substantially simplified derivations. Hence, this reformulation of the theoretical basis of the GRM is quite 97 useful for educational purposes. However, a practical benefit related to refractor depth pro file construction is also realized. The equations developed here allow point depth estimates of the interface to be calculated from the generalized time-depth function. Computation of a refractor depth profile then reduces to an interpolation problem (treated in Chapter 6). In contrast, Palmerβs time-depth function yields a circular locus of possible refractor positions. Depth profile calculation then involves the more complicated task of constructing an envelope to a set of circular arcs. As indicated above, the GRM is founded upon a two dimensional earth model with plane interfaces. A traveltime inversion technique that explicitly incorporates nonpiane refracting horizons into the model should yield more accurate results. Hence, the second portion of. this chapter proposes a 2D head wave traveltime interpretation method tentatively named critical offset refraction profiling. This inversion technique accommodates undulations in both the refracting interface and the surface, as well as horizontal variations in the velocity of the refracting layer. Detection and definition of such structural and velocity variations are subjects of current research (Palmer, 1991). Critical offset refraction profiling yields a mathematically exact solution to the 2D seismic refraction problem; point values of interface depth, interface dip, and refractor velocity are obtained. A continuous refractor depth profile can then constructed by interpolation. The dip estimates provide additional constraints on the interpolant. 5.2 Specialization to a 2D model Prior to deriving the GRM analysis functions, the head wave traveltime equation (3.23) must be specialized to a two dimensional earth model. In the three dimensional model, the orientation of interface i is described by a dip angle qfj (0 Γ§& <ir/2) and an azimuth angle S (0 S < 2ir). Note that the dip angle 4βj is non-negative. In terms of these two angles, the unit normal to the interface is n1 = (sin gj cos O)i + (sin Γ§t sin 02)j + (cos (5.la.) 98 and the vertical depth of the interface is z(x, y) = z(O, 0) β tan Γ§&j (x cos 8 + y sin 8). (5.lb) The vertical thickness of layer i (bounded from above and below by interfaces i and i + 1, respectively) is h(z, y) = h(0, 0) + x (cos O tan Γ§5j β cos8i tan qi+i) + y(sin Stan 4j β sin Sj.4 tan (5.lc) If the earth model is two dimensional, then all interface azimuth angles 6 are restricted to the two values 0 or ir. In this case, the above equations reduce to n = (sin Γ§Sj cos 8)i + (cos Γ§j)k, (5.2a) z(x) = z(0) β cos S tan 4j, (5.2b) h:(z) = h(O) + x(cos9tanq5 β cos tanc5i+1). (5.2c) Thus, the normal n1 has no y-component, and the depth z;(x) and thickness h1(z) are independent of the y-coordinate. For subsequent analysis of the two dimensional refraction problem, it is convenient to repararneterize the earth model in terms of interface dip angles that may be positive or negative. Hence, let the symbol Oj refer to an interface dip angle in a two dimensional layered earth. cpj is an acute angle (0 I < 7r/2) measured with respect to the +x axis; it is considered positive (negative) if the angle opens in the clockwise (counterclockwise) sense. This angle is related to the 3D interface orientation angles via = β (cos 9), where the azimuth 8 is restricted to the two values 0 or ir. Rewriting this relation as = βΓ§oj (cos 9i) and then substituting into equations (5.2a,b,c) gives the familiar expressions 99 = (β sin p)i + (cos y)k, (5.3a) z1(z) = zj(O) + x tanyj, (5.3b) h(x) = h(O) + sin(cpj l β (53) COS COS Γ§o The head wave traveltime formula relevant to GRM analysis is equation (3.23). This is specialized to a two dimensional model by setting ys = YR = 0, n = β 5111βpk, and k,z = cos cΒ°k. Equation (3.23) becomes T,(xs, XR) = β + β COS Yk + Pkk,z SilL Yk) (5.4) vkcoscok A further simplification is obtained by expressing the components of the unit propagation vectors Pik and q in terms of orientation angles. In the three dimensional situation, each of these vectors is described by two angles: Pik = (sin aik COS 7k)i + (sin aik Sβfl7jk)j + (cos a:k)k, (5.5a) qk = (sin/3 cos Sk)i + (Sfl/3jksin Sik)J + (cosflk)k. (5.5b) ak and /3j are poiar angles measured from the +z axis (0 ak, I3ik ir), and and 8ik are azimuthal angles measured from the +x axis (0 7ik, < 2ir). If the recording prolile is oriented perpendicular to the strike of the subsurface horizons, then afl propagation vectors are confined to the xz plane. The azimuth angles equal 0 or r, and expressions (5.5a,b) reduce to Pik = (+SiflUk)i+(C0Sak)k, q:k (+sin/3k)i+(cos/3k) . (5.6a,b) 100 The propagation direction of the critically refracted raypath segment is obtained by straight forward geometric analysis: Pkk = Β± [(cos cok)i + (sin Γ§ok)k], (5.6c) where the pius sign is used for R > xs and the negative sign for R < s. Finally, substituting equations (5.6a,b,c) into (5.4) gives the remarkably simple result h(xS)cosak β h(xR)cos/3k IZR β sI Tk(XS,XR) = + . (5.7) i=1 Vk COS Pk Equation (5.7) is a novel expression for surface-to-surface head wave traveltime that forms the basis for an alternative development of the GRM. It differs from the analogous equation given by Palmer (1980, p. 5) in several important ways: (i) layers are characterized by vertical thicknesses below source and receiver, (ii) raypath orientation angles are measured with respect to the vertical, (iii) the coefficient multiplying the source-receiver offset distance depends only on critical refractor quantities (i.e., velocity vk and dip and (iv) the earthβs surface may be nonhorizontai. These attributes facilitate a straightforward derivation of the various GRM analysis tools. The raypath angles in equation (5.7) depend on the recording proffle azimuth , which in turn is restricted to the two values 0 or r. However, raypath reciprocity requires that ak(7r) = βI3k(0) and f?(ir) β ak(0). These expressions imply that cosak(7L-) = βcos/3k(0), cosf3k(1r) = βcosak(0). (5.8a,b) Using these relations, it is easy to demonstrate that the traveltime formula (5.7) satisfies source-receiver reciprocity: Tk(XS,XR) = 101 5.3 Derivation of GRM parameters The two dimensional earth model used to derive the GRM analysis functions consists of uniform velocity layers bounded by plane, dipping interfaces. Figure 5.1 depicts the relevant critically refracted raypaths. Although the subsequent mathematical development assumes a multilayered earth, this single layer model is useful for illustrative purposes. Receivers are deployed between two sources Si and S2 with horizontal coordinates xs1 and Zs3, respectively; 1? denotes a typical receiver with coordinate ZR. The velocity analysis function and the generaiized time-depth function are constructed from the traveltimes of forward and reverse propagating head waves that arrive at the positions ZR + L/2 and ZR β L/2, respectively. L is an adjustable separation distance. If Β£ = 0, then arrivals at the single receiver at ZR are examined, and the GRM reduces to the conventional reciprocal method. GRM analysis seeks to determine the particular separation Β£ such that the upgoing segments of the two raypaths depart from the same point on the critically refracting horizon. This determination is made by evaluating the velocity analysis function and the time-depth function over a range of separation distances L, and then interpreting the plotted curves. The interpretive criteria that are used for selecting the appropriate distance Β£ are described by Palmer (1980, 1981, 1990, 1991) and are not investigated here. 5.3.1 Velocity analysis function The velocity analysis function for layer k is defined as T(xR;L) [zsix + L/2) β Tk(z52,xn β L/2) + Tk(xSizSa)]. (5.9) The last term on the right hand side is the reciprocal time (the shot-to-shot traveltime). 102 Fig. 5.1. Critically refracted raypaths used to derive the GRM analysis functions. R denotes a receiver located between the two surface sources S1 and S2. Substituting the head wave traveltime formula (5.7) into this definition yields kβi vi i=1 β h(:s2)[cosaik(ir)+cos/3ik(O)] (5.10) 0 XS1 L XR V2 + h(xj β L/2) COs/3jk(lr) β h(zR + L/2) COS/3k(0) + β xs1 2v t7k cos Γ§o i=1 where the dependence of the raypath angles on recording profile azimuth is explicit. Using equation (5.3c), the layer thicknesses at the horizontal coordinates XR Β± L/2 are written as h(xR Β± L/2) = h(xs1)+ (XR β Xs1 Β± L/2) sin(yj1 β cos Pi+1 cos cΒ° 103 Substituting this expression and relations (5.8a,b) into (5.10) and reducing yield the result T(xR;L) = -_(h(xs1)+ sin(j+1βcoj) cosak(0)βcos)6k(0) v 2 cos Γ§oj cos cpj j 2 +(ZRβxsjt 1 β sin(yi+iβy) cosak(0)+cosI3k(0)β1 (5.11) VCOSOJ v coscpjjcosoj 2 j - Hence, for a model consisting of uniform velocity layers bounded by plane interfaces, the velocity analysis function is a linear function of the distance (ZR β xs1). The slope of this line is used to estimate the velocity of the critically refracting layer. The intercept has a simple interpretation that is given in section 5.3.5. 5.3.2 Apparent horizontal velocity The apparent horizontal velocity of layer k is defined as 1 = 8T(z; L) (5 12) Va ΓZR Using expression (5.11) for the velocity analysis function immediately results in = 1 β i [sin(Yi+1_ci)] cosak(O)+cos31k ) (5.13) Va vj cos v1 cos Γ§o cos Oj 2 Thus, the apparent horizontal velocity of a refracting layer is influenced by all velocities and dips in the overburden, as well as by the dip of the critically refracting horizon. How ever, for dip angles commonly encountered in refraction exploration, formula (5.12) yields a reasonably accurate estimate of vk. A quantitative understanding of the accuracy of (5.12) can be obtained by analyzing a specific numerical example. Palmerβs (1980, p. 10) shallow subsurface model is quite useful for this purpose. Since the interfaces have unusually large 104 dips, this model provides a strong test of the utility of (5.12) for velocity estimation. Fur thermore, a comparison with the numerical values of Va computed via Palmerβs (1980, p. 9) formula serves as a check on the correctness of the current derivation. Palmer (1980) examines a four layer model defined by the parameters zi(0) = 0.00 m, = 0Β°, = 1000 m/s, Z2(0) 20.31 m, = 10Β°, V2 1800 mIs, Z3(0) = 59.55 m, O3 = 5Β°, = 2500 m/s, z4(0) = 80.25 m, cp = β15Β°, V4 = 4000 rn/s. (He does not give the vertical depths to the interfaces, but these may be derived from the specified geometric parameters). In order to evaluate equation (5.13), the raypath angles ak(O) and/3k(0) must be determined. Application of Snellβs law of refraction at interface i yields the two recursion relations βiβ1,k = sin1 [!z1 sin(ak + co)J β 13iβ1,lc = ir β sin [Vi_1 sin(/3k+ oi)j β Yi. (5.14) Starting with akk(0) + Yk =/3kk(O) + y = ir/2, these relations are evaluated successively fori=k,kβ1,...,3,2toobtain a2 = 43.7490Β°, /312 = 156.2510Β°, = 26.6536Β°, /13 = 160.9118Β°, a23 = 41.0545Β°, /323 = 128.9455Β°, a1 = 22.5156Β°, /314 = 174.7119Β°, a24 = 32.9588Β°, /324 = 161.6659Β°, = 53.6822Β°, /334 = 156.3178Β°, Finally, using these numerical values for the raypath angles in formula (5.13) yields k=2: Va1827.77m/s, k=3: Va=2578.59m/s, k=4: Va4209.65m/s. 105 These calculated apparent velocities agree closely with Palmerβs (1980, p. 11) values, and thus confirm the validity of equation (5.13). Moreover, it is evident that Va is a good ap proximation to the true refractor velocity v; percentage errors are reasonably small (+1.5%, +3.1%, and +5.2% for layers 2, 3, and 4, respectively). As Palmer (1980, p. 9) points out, an improved estimate of refractor velocity can be obtained if the dip Γ§o of the critically refract ing interface is known. This estimate is derived as follows. Assume that the dip angles of all interfaces overlying the critically refracting horizon vanish (pj = 0 for i = 1,2,. . . , k β 1). Then, equation (5.13) becomes 1 1 β tan yj [cos akβ1,k(O) + COS /3kβi ,k(0) Va vkcoscpk Vk_i L 2 Geometric analysis reveals that the raypath angles within the layer immediately above interface k are given by ak_i,k(0) = uk β iI and13k_1,k(0) β uk + cpkl, where = sin(vk_i/vk) is the critical angle at the kth interface. Substituting these relations into the above expression for Va and reducing yield Va vk (5.15) Cpjβ Thus, under the stated assumptions, the apparent velocity exceeds the true refractor velocity. Equation (5.15) holds in an approximate sense if the dip angles of all interfaces overlying the critical horizon are small, rather than identically zero. A dip corrected apparent horizontal velocity can be defined as V cos Pk Using the above model parameters, the following values are obtained: k=2: V=1800.00m/s, k=3: V=2568.77m/s, k=4: 1β=4066.21m/s. Obviously, Tβ is an improved estimate of the true refractor velocity; percentage errors are 0%, +2.8%, and +1.7% for layers 2, 3, and 4, respectively. Note that VL calculated for layer 2 is exactly correct, because the dip of the overlying interface (the surface) is zero. Interestingly, the error does not increase monotonically with layer depth. 106 As indicated previously, the condition for exact validity of equation (5.15) is the dip angles of all interfaces above the critical horizon vanish. Palmer (1980, p. 9) states that (5.15) is approximately correct if the differences in the dip angles between successive in terfaces (coj+1 β ep,) are negligible. However, this assertion is easily disproved. Expression (5.13) for the apparent horizontal velocity is evaluated for a multilayered model with veloc ities v = 1000 + 100 i rn/s and dips Γ§o iΒ°, ( i = 1,2,.. . ,50). All differential dip angles equal 10 and thus are small. A comparison between the apparent velocities and the true layer velocity is illustrated in Figure 5.2. In particular, both Va and V depart appreciably irom the correct velocity value as the layer index increases, i.e., as the dips of the bounding interfaces become large. This behaviour confirms that the approximation (5.15) requires small absolute dip angles, rather than small differential dip angles. Fig. 5.2. Comparison of apparent velocity Va and dip-corrected apparent velocity V with the true layer velocity vk, for a multilayered earth model. Another model which can be used to establish the veracity of (5.13) consists of a set of parallel, but dipping layers. Assume that all interfaces (including surface and critical refractor) have the same dip angle (Γ§o = Γ§o, for i = 1,2,.. . , k). Obviously, the differential U) 10 20 30 40 50 layer number 107 dip angles are identically zero. Expression (5.13) immediately yields the expected result: V0, = vk cos cp. In this case, the apparent velocity is less than the true refractor velocity. Finally, an alternate representation of the apparent horizontal velocity is derived from the slope/intercept formula for head wave traveltime: Tk(xs,X,P) = mk(J!)X + bk(xs). (5.16) Explicit expressions for the slope mk(βIβ) and intercept bk(xs) can be obtained by specializing equation (3.12) to a two dimensional earth model. Using this equation in the definition of the velocity analysis function yields 1 β 8T β mk(O)+mk(lr) 517 VaOR 2 β’ Thus, the apparent horizontal slowness is merely the arithmetic mean of the slopes of the forward and reverse traveltime curves. This implies that Va is the harmonic mean of the velocities measured on the forward and reverse spreads. Defining vf 1/m(O) and Vr 1/mk(Tr), then 1 11 1 VfV,. β = β β + β , or V0 = 2V0 2 Vf Vr Vf+Vr Lankston (1990) and Palmer (1990) state that the harmonic mean of v1 and v7 equals the true layer velocity. The present analysis reveals that claim to be erroneous. Rather, the harmonic mean V0 is an approximation to the true velocity vk. 108 5.3.3 Generalized time-depth The generalized time-depth function for interface k is defined as Td(XR; L) [Tk(zs ,ZR + Β£/2) + Tk(ZS2,ZR β L/2) β , xs2) β ]. (5.18) Substituting from equation (5.7) for the head wave traveltimes gives Td(ZR;L) [cosak(0)_cosik(0)] I, β1L j 1 Γ§- 1 s1n(+1 β yi) cosak(0) + cos8k(O) 1 + 2jvlccoScok 4vi CoSyjcoScpj 2 Va However, from equation (5.13) for the apparent horizontal velocity, the term in braces is identically zero. The above expression reduces to Td(ZR;L) = h(x1) [cosaik;cosf3ilC] (5.19) where the explicit dependence of the ray angles on recording proffle azimuth has been sup pressed because, from relations (5.8a,b), the difference cos k(I1)_cos/3(1P) is independent of azimuth. The time-depth of a critically refracting interface at coordinate ZR is related to the vertical thicknesses of all overlying layers at the same position. Plotted time-depth sec tions give a qualitative indication of refractor structure, in a manner analogous to reflection time sections. Conversion of the time-depth to actual depth requires velocity information, and is briefly discussed in the next section. 109 The slope/intercept formula for head wave traveltime also yields an alternative expres sion for the time-depth function. Substituting equation (5.16) into the definition of the generalized time-depth and reducing gives Td(XR;L) = [mk(lr) ; mk(O)] ( β XR) + bk(ZS2) The reciprocal time (or shot-to-shot traveltime) is Tk(zS1,as2) = Tk(as2,cs1). Thus, from - the slope/intercept formula mk(0)(x52 β xs1) + bk(x51) = mk(lr)(x52β z51) + bk(x52). Hence bk(zS1) β bk(xS2) mk(7r)βmk(O) = xS2 β Finally, substituting this result into the above expression for the time-depth function yields Td(cR;L) = bk(z51) Zs2 XR . bk( a52) XRS1 (5.20) 2 x2 β 2 Xs2 β Xs1 The time-depth function at receiver coordinate ZR is a linear interpolation of the half- intercept times at the two bounding source coordinates. Thus, for the chosen earth model, the GRM time-depth section is identical to the half-intercept time section used in classical refraction traveltime analysis (Wyrobek, 1956). 110 5.3. Depth conversion factor Expression (5.19) for the generalized time-depth may be rewritten as Td(XR;L) = h(x) (5.21) Vjk where the depth conversion factor is defined as Vji, 2v (5.22a) COS ajk β COSI3ik The depth conversion factoi has physical dimension of velocity. Since the raypath angles exile and /3ik are usually not determined in a GRM interpretation, Vk cannot be evaluated numer ically. However, various approximations to the actual depth conversion factor are useful. If the cosines of the raypath angles are replaced with the expressions appropriate for horizontal layering (cos exile β cos = β (Vj/Vk)2), then equation (5.22a) is approximated by VjVk = ___ ___ _ . (5.22b) Alternately, the true layer velocities in (5.22b) can be replaced by apparent horizontal ve locities to obtain - ai ale = _ _ , (5.22c) i/VVa where Vaj denotes the apparent horizontal velocity of layer i. Evaluating these expressions with the parameters defining the shallow four-layer model yields the following numerical values (in m/s): 111 i k V Vik 1 2 1221.23 1202.68 1194.66 1 3 1087.69 1091.09 1084.91 2 3 2603.66 2593.76 2591.17 1 4 1041.93 1032.80 1029.47 2 4 2013.08 2015.61 2029.00 3 4 3315.54 3202.56 3262.22 The approximations appear to be reasonable; the largest percentage error is β3.4% (for v34). The depth conversion factors are used to calculate layer thicknesses from the observed time-depths via expression (5.21). For k = 2, this equation yields hl(XR) = v12 Td(XR; L), (5.23a) while for k > 2: = vk_1,k [T4XR; L) β (5.23b) Vertical depths to the refracting interfaces are obtained simply by summing the layer thick nesses: zkfrR) = zl(XR) + hj(XR) for k 2. Note that the depth zk(XR) is an unambiguously defined function of the receiver coordinate ZR. In contrast, the βdepthβ cal culated via Palmerβs (1980) time-depth function is a radius of a circular arc centered at the receiver position. Continuing with the illustrative example, numerical values of the time-depth for the three subsurface interfaces can be obtained by evaluating expression (5.20). Sources are located at the horizontal coordinates zs1 = β100 m and xs2 = +60 m, and the receiver is positioned at the coordinate origin. The results are k = 2: Td(0;L) = 16.63ms, k = 3: Td(0;L) = 33.74ms, k = 4: Td(0;L) = 45.23ms. 112 Finally, vaiues of interface depths at ZR = 0 are calculated using the exact and approximate depth conversion factors. The results are (in m): z(Q) (0) 2 20.31 20.00 (β1.5%) 19.87 (β2.2%) 3 59.55 59.98 (+0.7%) 59.85 (+0.5%) 4 80.25 79.29 (β1.2%) 80.15 (β0.1%) The associated percentage errors indicate that, for this example, an adequate result is ob tained by using the approximate depth conversion factors. 5.3.5 Time-depth near a shotpoint The intercept on the time axis of the velocity analysis function (5.11) has a simple, but nevertheless useful, interpretation. From formula (5.3c), the vertical thickness of the th layer at the coordinate xs1 + Β£/2 can be written as hj(z1 + L/2) = h(xs1)+ sin(co:+1 β2 cos Pi+1 cos Γ§o Substituting this result into expression (5.11) for the velocity analysis function yields T1,(XR;L) β h(xs1+ L/2) [cosoik β COS/3jk] + ZR βZS1 However, from equation (5.19), the summation term is recognized as the generalized time depth at the coordinate xs1 + L/2. Thus T(xR; L) = Td(ZS1 + L/2; L) + Z R ZS1 (5.24) 113 The time axis intercept of the linear velocity analysis function equals the interface time- depth at a point offset a distance L/2 from source 1. This interesting relationship provides a tool for estimating refractor depths near a source, where head waves are either nonexistent or difficult to pick accurately. 5.4 Critical offset refraction profiling 5.4.1 Earth model The two dimensional earth model considered in this section is ifiustrated in Figure 5.3. A single layer with uniform velocity v is underlain by a medium where velocity v2(z) depends solely on the horizontal coordinate. Surface and subsurface horizons are undulating and are defined by the depth functions zl(x) and z2(x), respectively. The dip angles of these two interfaces are given by β1 Fdzi(z) dz2(x) = tan dx S02(Z) = tan dx These relations define the sign convention associated with interface dip. Forward and reverse traveltime data axe recorded by receivers distributed along the surface. The observed traveltime curves are designated Tf(z) and Tr(x), respectively. The slopes of the two traveltime curves are dTf(z) β dTr() mf(x) dx m(x) = β dx Note that a negative sign is included in the definition ofm7(x). From the measured slopes of the traveltime curves, the angles of incidence of rays arriving at the surface can be calculated. Let the functions (x) and ,t(x) refer to the incidence angles of the forward proffle and reverse profile raypaths, respectively. Then, straightforward geometric analysis yields the expressions 114 Fig. 5.3. A 2D earth model with undulating surface and refractor topography. Forward and reverse raypa.ths exiting from the same point C on the critical horizon intersect the surface at B and A, respectively. M is the midpoint coordinate of A and B. ILf(Z) = sink cosS1(x)mf(a)] + yi(x), = sin [vlcoscol(x)mr(x)] β jtffr) and are acute angles measured with respect to the vertical. The previously established sign conventions for surface dip and traveltime slope determine the positive sense of these two incidence angles. In particular, jLf(c) is positive if the angle opens in the clockwise sense. The converse holds for ILr(x); it is positive if the angle opens in the counterclockwise sense (see Figure 5.4). 0 Lc(xM) XM zi A B vi V2(X) 115 The measured traveltime data pertain to waves that are refracted at the subsurface horizon z2(r). In general, the raypaths of these waves may undergo noncritical refractions at the interface. However, the inversion procedure discussed here adopts the following fun damental assumption: the critically refracted portion of the total raypath is coincident with the undulating subsurface interface. Although this is an approximation, it is sufficiently accurate if the structure on the interface is not too severe. In the presence of syncinal structure, the approximation is reasonable because there is a tendency for the diffracted raypath to follow the interface (see, for example, Figure 7.1). Since raypaths penetrate be -neath anticinal structures, the approximation is less accurate in these cases. Nevertheless, refraction seismologists commonly adopt this assumption as a useful working hypothesis. In fact, it underlies the popular interpretive rule known as the βlaw of parallelismβ of refraction traveltime curves (Rockwell, 1967; SjΓΆgren, 1980; Ackermann et al., 1986; Brflckl, 1987). Figure 5.3 depicts forward and reverse raypaths exiting from the same point C on the subsurface interface. If i2(z) sin1 [vl/v2(x)] is a variable critical angle defined along interface 2, then these rays depart at an anglei2(xc) with respect to the normal to z2(2). Points of intersection of the reverse and forward rays with the surface zi (x) are designated A and B, respectively. These two points are separated by the critical offset distance, i.e., the minimum horizontal distance between a surface source and receiver at which a critical refraction can be observed. Obviously, the critical offset varies along the line of proffie due to the changing depth and dip of the two interfaces, as well as the the varying velocity v2() of the lower medium. If XM denotes midpoint coordinate of A and B, then a critical offset function L(xM) can be defined at every midpoint position on the surface. The inversion technique outlined below is based on the fact that an exact solution to the refraction problem is possible if this function is known. A procedure for estimating the critical offset function from the measured traveltime data is discussed in a subsequent section. 5.4.2 Inversion method If the two points A and B on the surface are separated by the critical offset distance, then a closed form solution for the subsurface model parameters can be obtained. Figure 116 5.4 indicates that β XBβXC tan ttf(XB) β zC β Z5 XC β X4 tan1t(xA) β zC β Z Solving these equations for the two unknown coordinates (XC, zc) yields XC = zC = (XB β ZA) COS /Lf(XB) cos liT sin (ti.e() + (5.25a) (5.251,) These relations establish the position of the point of critical refraction C in terms of known quantities. Furthermore, since the sum of the angles of triangle ABC equals ir radians, the critical angle at C is given by pf(XB) + lLr(XA) i2(XC) = 2 (5.26) Hence, the refractor velocity at point C is determined by v2(xC) = vi! sini2(XC). Finally, the local dip of the refracting interface is calculated by analysis of either triangle ACD or triangle BCE in Figure 5.4. The result is S02(XC) = lif(XB) β XB COS lif(XB) sin iir(ZA) + XA Sifl ILf(XB) cos Pr(XA) sin (lif(XB) + (β β zA)sinl&f(XB)sinlir(ZA) sin (lLf(XB) + 1Lr(2A)) ZB sin pf(XB)cosILD(ZA) + zAcoslif(XB)sinlLr(XA) sin (lif(XB) + ltr(XA)) 2 (5.27) 117 Fig. 5.4. An expanded view of triangle ABC in Fig. 5.3. Angles j and p describe the orientation of raypath segments CB and UA relative to vertical, respectively. The critical and dip angles at point C on the refractor are i2 and Y2, respectively. Expressions (5.25) through (5.27) constitute a mathematically exact solution of the specified seismic refraction problem. The formulae are not restricted to the particular geometric situation depicted in Figure 5.4, but remain valid if the refracting point C does not reside between the surface points A and B (i.e., if either pf(zB) <0 or ir(XA) < 0). The above inversion formulae for critical angle and dip are similar to the classical re fraction solution for an earth model with plane interfaces and uniform velocities. In that case, the forward and reverse traveltime curves are straight lines with slopes that are inde pendent of the horizontal coordinate x. Hence, the forward and reverse angles of incidence are also constant: ILf(x) = and p7(x) = pr. Equations (5.26) and (5.27) reduce to the familiar expressions i2 = (ELf + p1)/2 and Y2 = (jLf β ,u)/2. There are no analogues to A 118 equations (5.25) in the classical solution, because vertical depths to the refracting interface are determined only at the source locations. A set of point estimates of the position of the critically refracting horizon is obtained by repeated application of formulae (5.25a,b) at various midpoint locations along the line. A profile of the interface can then be constructed by interpolating these points. Additional constraints on. the interpolant are provided by the interface dip estimates at each refract ing point C. Assuming an error-free inversion, the constructed depth proffle should pass through the coordinates (ac, zc) and have slope tan cp2(xc) at this point. A procedure for constructing a smooth interpolation of a set of point estimates of depth and dip is presented in Chapter 6. 5.4. 3. Critical offset determination The inversion method outlined above requires the critical offset distance L(xM) to be known at each analysis midpoint along the profile line. Jones and Jovanovich (1985) describe a computational technique for estimating critical offsets at the source positions on a line. This section discusses an extension of their method to arbitrary locations along a reversed refraction proffle. Moreover, an analytical, rather than computational solution to the problem is given. Consider two surface locations with a common midpoint coordinate XM and separated by the horizontal distance L. A candidate critical refraction raypath can be constructed from these two points by projecting the incoming raypath segments downward until they intersect (Figure 5.5). The projection directions are determined from the known incidence angles f(rM + L/2) and (ZM β L/2). The propagation time along this hypothetical raypath is called a βpredictedβ traveltime and can be derived via geometric analysis. It is 119 T7,,.d(ZM;L) =. Lc(xM) COS1Uf(XM + L/2) + co5r(xM Β£12) (Itf(XM + L/2) + Itr(XM β Β£12)) zl(ZM + Β£72) β zl(XM β Β£12) Sin1Uf(ZM + L/2) β sin/ir(aM β L/2) + L/2) + i?(xM β Β£12)) Note that this expression simplifies considerably if the surface is horizontal (i.e., zi(x) = zi(O) = constant). At the same two surface locations, an βobservedβ traveltime is determined from the measured forward and reverse traveltime curves as follows: XM A \1/ B Z2(X) Fig. 5.5. Determination of the critical offset distance L(xM) at the midpoint coordinate ZM. Raypath segments projected downward from A and B on the surface intersect at point C on the refracting interface. T0b3(XM; L) = Tf(XM + L/2) + Tr(ZM β L/2) β TR, (5.29) 120 where TR is the reciprocal time: TR = Tf(xS2)= Tr(xj). For a fixed midpoint coordinate XM, expressions (5.28) and (5.29) are evaluated over a range of values of the separation distance L, extending from well below to well above the expected critical offset distance. Figure 5.4 indicates that whenT7,.d(xM; L) =T0b3(XM; L), then L = L(aM). The common time calculated at this separation distance is the critical offset traveltime Tc(XM) associated with the midpoint position. Application of this analysis technique at various midpoint locations along the profile line yields several point estimates of the critical offset function L(x). 5.5 Conclusion The theoretical foundation of the generalized reciprocal method of refraction traveltime inversion has been reformulated in terms of layers characterized by vertical, rather than normal, thicknesses. This reformulation results in a simplification of the derivations of the two GRM analysis functions, and thus has educational value. However, it also allows point estimates of interface depth to be calculated from the measured traveltimes. A refractor depth profile can then be calculated by interpolating the depth points. Although the GRM is based on a simple 2D model with plane interfaces, interpreters apply the technique in more realistic environments with nonpiane refracting horizons. The results that are obtained should then be viewed as approximate. In contrast, the proposed method of critical offset refraction profiling explicitly incor porates undulating interfaces and variable refractor velocities into the model. A mathe matically exact solution to the traveltime inverse problem indicates that point estimates of interface depth, interface dip, and refractor velocity can all be obtained from the observed arrival times. These point estimates are positioned properly in space; there is no need for a subsequent migration of the solution. In common with many other refraction interpretation methods, independent knowledge of the overburden velocity v is required. Although the derivations in this chapter assume a uniform layer velocity, this is not a practical necessity. Rather, it is only necessary that vi be reasonably constant over the extent of the triangle ABC in Figure 5.3. Lateral variations with a characteristic scale greater than the horizontal i2 1 dimension of the triangle are allowed. The primary approximation adopted by this technique relates to the refracted raypath: it remains coincident with the undulating interface. The implications of this approximation for the inversion should be investigated more thoroughly. Finally, extension of the method to multilayered earth models represents another avenue of developement. Using the present algorithm, this could be achieved via a layer stripping approach. 122 CHAPTER 6 CONSTRUCTION OF A SMOOTH REFRACTOR DEPTH PROFILE 6.1 Introduction Linear inverse theory provides techniques for constructing acceptable models that are consistent with the available data. However, if acceptabffity is judged by the data misfit alone, then there are usually infinitely many equally valid solutions. A specific solution can be obtained only by extremizing some functional of the model, subject to a requirement that the final model is in satisfactory agreement with the observations. Since different functionals can produce models with significantly different character, an essential step in the solution of any inverse problem is to select the βrightβ functional. Minimization of the 12 norm of the model yields a βsmallestβ model that often has high variability and hence may be difficult to interpret. Similarly, models can be constructed by minimizing the 12 norm of the first or second derivative; these are often referred to as βflattestβ or βsmoothestβ models, respectively. These latter types are particularly useful because they can be interpreted as minimum structure solutions for many geophysical inverse problems. It is anticipated that the true earth model that gives rise to the measured data has at least as much structure as observed on these constructed models. Methods for calculating minimum structure solutions have been described by various authors. These methods differ in the manner in which the required additional information regarding the model is treated. Johnson and Gilbert (1972) construct a smoothest model by explicitly incorporating endpoint values of the model and its derivative into the objective functional to be minimized. Oldenburg (1984), in computing a flattest model, specifies an a priori model value at an endpoint. Parker, Shure, and ilildebrand (1987) minimize a seminorm of the model. This chapter discusses two generalizations to the conventional method for constructing a smoothest model. In this technique, the governing equations are integrated by parts twice; each such integration introduces extra parameters into the model construction problem. Numerical values for these additional parameters are obtained either by (i) prescribing a priori values at any abscissae within the interval of model definition, 123 or (ii) optimizing the appropriate objective function with respect to the parameters. These generalizations are readily extended to model construction problems that involve minimizing the 12 norm of other derivatives of the model. In this chapter, the smoothest model construction formalism is applied to the problem of determining a depth proffle for a critically refracting horizon. It has only recently been demonstrated that refraction traveltime data can be processed to yield point estimates of interface depth (Jones and Jovanovich, 1986; Pavienkin et aL, 1986; All Ak, 1990). Moreover, as shown in the previous chapter, local dip estimates can also be derived from the traveltimes. -A refractor depth proffle can then be calculated by smoothly interpolating all of these data. This procedure provides a useful alternative βto the popular time-to-depth conversion method associated with the GRM (Hatherly, 1980; Hatherly and Neville, 1986; Kilty et al., 1986; Lankston and Lankston, 1986; Lankston, 1989; Palmer, 1990, 1991). In the GRM technique, circular arcs -with radii given by expression (5.23a) are ceutered at each analysis location along the seismic line. The position of the refracting interface is then estimated by the envelope of the set of arcs. This envelope is approximated by tangent lines to pairs of adjacent arcs; problems specific to this method are identified and discussed by Hatherly (1980). The technique has antecedents in the depth determination procedures of Thornburgh (1930), Dix (1941), Tarraut (1956), Hales (1958), and Hawkins (1961) where possible refractor loci, rather then point depth estimates, are calculated from the traveltime data. The interpolation method described in this chapter circumvents those problems related to calculating an envelope of a set of circular arcs. It also possesses several specific advantages regarding treatment of the data. However, perhaps the most important difference pertains directly to the derived models: the method proposed here yields the smoothest refractor depth profile consistent with the given data. This attribute is not shared by the model constructed via the GRM technique. A general theoretical developement of the smoothest model construction method is pre sented before discussing the particular application to interpolation. Hence, the results in this chapter are useful for the solution of other geophysical inverse problems that may require a 124 smooth model. Assume that an N-vector of observed data e0b3 is related to the true earth model mj(x) via a Fredholni integral equation of the first kind: b e0b3 = j g(x) mt(x) dx + Se. (6.1) g(z) is a vector of known kernel functions and Se is a vector of additive random errors. The true earth model mt(x) is unknown and is to be estimated. Predicted data generated by a constructed earth model m(x) are given by the linear functional b ed(m) j g(x) m(x) dx. (6.2) In the following development, the 12 norm is used to quantitatively measure both model structure and data misfit. For example, the square of the 12 norm of the second derivative of the model is b )j rnβ 112 = j mβ(x)2 dx. Similarly, the square of the 12 norm of the misfit between observed and predicted data is given by II e03 β e,,.d(m) 12 = [eobs βe1.d(m)j [cobs β elfl.d(m)}, where the superscript T indicates transposition. 6.2 Smoothest model construction A detailed discussion of the method for calculating the flattest model is given by Aldridge et al. (1991). The analogous derivation for the smoothest model is conceptually similar, although algebraically more complicated. 125 6.2.1 Modified data equation Initially, equation (6.2) is integrated by parts to obtain a modified data equation satisfied by the first derivative of the constructed model: b ed(rn) m(b)h(b) β j h(x)mβ(x) dx, (6.3) where h(x) is the indefinite integral of the original kernel function vector g(x): h(x) = / g(u) du. (6.4) The value of the model at the endpoint b is required to evaluate the right-hand side of equation (6.3). However, the expression is easily rewritten in terms of a model value m(ci) specified at an arbitrary abscissa c within the closed interval [a, b]. From the fundamental theorem of calculus b b m(b) β m(ci) = j mβ(x) dx = j H(x β Cl) mβ(x) dx, (6.5) where H(x) is the Heaviside unit step function. Eliminating m(b) from equations (6.3) and (6.5) yields ed(m) = m(ci) h(b) β / {h(x) β h(b)H(x β ci)] mβ(x) dx. (6.6) In the model construction problem, m(ci) is considered to be a parameter independent of the abscissa Cl. This independence is emphasized by writing d1 for m(ci). Equation (6.6) becomes b, ed(mβ) = d1h(b) β j p(x;cl)mβ(x) dx, (6.7) 126 where a new kernel function vector is defined by p(x; ci) h(x) β h(b) H(x β ci). (6.8) Note that the predicted data are now considered to be a functional of the first derivative mβ, rather than the model m itself. Furthermore, due to the presence of the term d1h(b), these data are a nonlinear functional of mβ. If the above procedure is repeated, then a modified data equation satisfied by the second derivative of the constructed model can be derived. Thus, expression (6.7) is integrated by parts and the fundamental theorem is used to transfer endpoint information to another arbitrary location C2 with the interval [a, bJ. The result is ed(mβ) = d1 h(b) β d2 [k(b) β (b β ci)h(b)j + J p(x; ci, c2) mβ(x) dx, (6.9) where the kernel function vector is defined by p(x;cl,c2) k(z) β h(b)R(z β ci) β {k(b) β (bβ ci)h(b)JH(x β c2). (6.10) ci and c2 are two abscissae in [a, b} where the parameters d1 m(ci) and d2 mβ(c2) are defined. Also, k(x) is the indefinite integral of h(z): k(x) = / h(u) du. (6.11) Finally, R(x) in equation (6.10) is the ramp function with unit slope: R(z) = xH(x). Expression (6.9) is the desired result. It indicates that the predicted data are a nonlinear functional of the second derivative of the constructed model. In the particular case where both c1 and C2 are restricted to the upper endpoint b, formulae (6.9) and (6.10) reduce to familiar forms. 127 6.2..2 Objective function In order to construct a smooth model that simultaneously minimizes the misfit between observed and predicted data, consider the objective functional mmβ 112 + βN [eobi β eWd(mβ)] 112 . (6.12) The scalar jt (0 < +oo) is a tradeoff parameter that controls the relative importance of the two terms. Equation (6.12) also includes a model derivative weighting function w(z) and a data weighting matrix W. The function w(z) can be used to emphasize certain portions of mβ(z) during subsequent minimization of 4β. Similarly, the matrix W aiiows a point-by- point weighting of the input data according to prescribed criteria. Commonly, W2 is taken to be the inverse of the covariance matrix of the observational error 6e. However, the matrix in (6.12) is not restricted to this specific type of weighting. Appendix D demonstrates that the second derivative function that extremizes 4 must be a linear combination of the weighted kernels pj(x; c, c2)/w(x) j = 1,2,. . . , N: mβ(z) = aT p(z;ci,c2) (6.13) where a is an N-vector of coefficients. Using this representation for mβ(z), it is straightfor ward to demonstrate that wmβ 112 = aTF(ci,c)a (6.14) and e,,.d = d1h(b) β d2 [k(b) β (bβ ci)h(b)] + r(ci,c2)a, (6.15) where F(ci, c2) is an inner product matrix formed from the kernel function vector p(x; c1, c2): b F(c1,c2) j w(x)2p(x;cl,c2) p(x;ci,c2)T dx. (6.16) 128 The matrix F(c1,c2) is symmetric and, if the original kernel functions gj(x) are linearly independent on [a, bj, positive deilnite. Substituting expressions (6.14) and (6.15) into (6.12) converts the objective functional into a quadratic form in a, d1, and d2: (a,ci,c2,cli,d)= aTr(ci,c2){WTWr(ci,c+,ttI]a + 2 [d1h(b) β d2 [k(b) β (b β c )h(b)] β eobs]TT r(ci, C2) a + dh(b)TWTWh(b) β 2dldh(b)TWTW[k(b) β (bβ ci)h(b)] + 4 [i β (b β cl)h(b)]TWTW [k(b) β (b β ci)h(b)] β 2d1 h(b)TVSITWe0b3 + 2d [k(b)_(b_ci)h(b)jTv.rTwe0,3 +e3WTWe0b3. (6.17) Dependence on the two abscissae ci and cz is via the inner product matrix Iβ(ci, c2) (as well as explicit dependence on Cl). Note that expression (6.17) is a complete quadratic form in the independent variables a, d1, and d2 in the sense that afl possible terms are present. In contrast, the analogous quadratic objective function given by Johnson and Gilbert (1972) lacks the cross product term in d1 and d2, the linear terms in both d1 and d2, and the constant term. Additionally, the remaining terms differ in detail from those given above. 6.2.3 Extremizing the objective function Straightforward methods of multiva.riable calculus are now used to extrernize the objec tive function with respect to all of the variables relevant to the problem. Thus, calculating O/8a and setting the result equal to 0 yields the linear system {F(ci, c2) +11Γ§VVTW)_l]a = e0 β d1 h(b) + d2 {k(b) β (b β ci)h(b)]. (6.18) 129 If values for the four scalars (Cl, d1, C2, (12) are prescribed, then (6.18) can be solved for the coefficient vector a. This offers an interesting alternative to the conventional method for constructing a smoothest model, where Cl and c2 are typically restricted to one of the end points a or b. However, for the current analysis, optimum values of the four parameters are desired. Thus, extremizing β with respect to d1 and d2 yields the two additional expressions [wTWh(b)] [di h(b) β d2 {k(b) β (b β ci)h(b)] β e0b + F(c1c2)a] = 0, [WTW [k(b) β (b β cl)h(b)J] [dlh(b) β d2 [k(b) β (b β ci)h(b)] β e0b3 + r(c1,c2)a] = 0. Substituting from equation (6.18) immediately reduces these expressions to the simpler forms aTh(b) = 0, (6.19) aTk(b) = 0. (6.20) Geometricafly, these conditions imply that the N-dimensional coefficient vector a must be orthogonal to the 2-dimensional subspace spanned by h(b) and k(b). Finally, the objective function 4β must be extremized with respect to the two abscissae Cl and C2. However, Appendix E demonstrates the remarkable result that the derivatives Γ΄4β/Ocl and Γ΄4β/8c2 both vanish if the prior conditions (6.18), (6.19), and (6.20) hold. Thus, these three equations are sufficient to extremize 4β with respect to the live quantities (a , Cl, d, C2, d2). As long as these conditions are imposed, any convenient pair of abscissae Cl and C2 can be chosen to evaluate the inner product matrix Tβ(Cl, c2). Expressions (6.18) through (6.20) are then a system of N + 2 equations in the N + 2 unknowns a3, d1, and d2, and can be solved simultaneously with standard techniques of linear algebra. However, 130 additional insight is obtained by first eliminating the coefficient vector a to form a simple 2 x 2 system for d1 and d2: Ad = b. (6.21) In this expression, d = (d1,d2)T, and the elements of the matrix A and the vector b are given by all = h(b)Tr_(jt) h(b), a = (b β cl)all β k(b)TF_l() h(b), a2l = h(b)Tr_l(t) k(b), a22 = (b β cl)a21 β k(b)Tr(t) k(b), = h(b)Tr(I)eob$, = k(b)TF_l(,I)eobs. (6.22) Here rβ1(,t) stands for {r(ci, c2) +,t(WTW)_hj . If the special cases h(b) = 0, k(b) = 0, and k(b) = /3 h(b) (where /3 is a constant) are excluded, then it can be shown via the Cauchy Schwartz inequality that the determinant of this system is nonzero. The optimum values of d1 and d2 obtained by solving (6.21) are linear combinations of the observed data eobs. After these two parameters are determined, the coefficient vector a is found by solving equation (6.18). 6.2.4 Constructing the model The model m(x) is obtained by integrating the second derivative mβ() in equation (6.13) twice. Integration constants are chosen so that the predicted datae3,,.d(m) reproduce the observed data e0b3 in the limit as the tradeoff parameter JL approaches zero. Thus m() = di+d2(x_ci)+aTj jp(u;cic2)du dv. (6.23) Optimizing the values of the parameters d1 and d2 in the above manner has beneficial implications for the curvature of the model. In general, the second derivative of the model 131 is given by expression (6.13). Substituting in the explicit form for the kernel function vector p(x;cl,c2) yields mβ(x) = w(x)_2{aTk(x) β aTh(b) {R(x β ci) β (b β ci)H(z β c2)j β aTk(b)H(x β c2)}. If conditions (6.19) and (6.20) hold, then the second derivative simplifies considerably to mβ(z) = a Tk(z)/w(x)2 These extra conditions remove a step discontinuity in the second derivative at a = cz due to the Heaviside function. Assuming that k(x) and w(z) are continuous, then the constructed model is twice continuously differentiable, i.e. m(x) E C2 on (a,b). This model is also unique except in those cases where the determinant of the 2 x 2 matrix A vanishes. Three particular cases are readily identified by analysis of det A. From equations (6.22): det A = [hbT rββ() k(b)] 2 β [hbT T() h(b)j [k(b)Tr_1(tt) k(b)j. Obviously, this determinant vanishes if either h(b) or k(b) equals 0. The sole remaining case arises in the following manner. Since () is positive definite, rββ/2(p) exists. Define two vectors u and v as: rβ/2()k(b). The above expression for the determinant becomes det A = (uT v)2 β (uT u) (vT v). If u /3v (where /3 is a constant), then the Cauchy-Schwartz inequality implies that u 112 v > (uTv)2. Hence, the determinant is nonvanishing in this situation, which is equivalent to the condition h(b) /3k(b). In these three cases, the constructed model is unique only to within an arbitrary additive constant ao or linear function a + aix. These situations are now described in detail. 132 Case 1: h(b) = 0. This situation arises if the original kernel functions g3(x) all possess zero area (e.g., Oldenburg, 1981). The system Ad = b reduces to a single equation that can be solved for d2: L-(ITpβ1( d β e03 2 β β k(b)Iβ-1L)k(b) Also, the linear system for the coefficient vector a simplifies to [r(ci, C2) + (ITT)IJ)_1ja = e0,3 + d2 k(b). Since the right hand side is known, solution for a is possible. The model m(x) is then calculated from (6.23) by picking any value for the parameter d1. Thus, constructed models differ by an arbitrary additive constant. Case 2: k(b)= 0. The 2 x 2 system Ad b reduces to a single relation between d1 and d2: h(h\TTββ1( .. e0b3 β 1 + β Cj 2 = h(b)TT ()h(b) = Also, expression (6.18) becomes [r(ci,c2) + ifVVTW)hja = e0,3 β It is again possible to solve for the coefficient vector a. The model is constructed via (6.23) after picking a value for either d1 or d; the remaining parameter is obtained from the above relation between d1 and d2. Constructed models then differ by an arbitrary linear function of x. Case 3: k(b) = /3 h(b). This case is a generalization of the previous one. The two equations in (6.21) are no longer independent, but reduce to the single expression d1 + (bβci β/3)d2 = 133 Once again, equation (6.18) becomes [F(cl,c2) + t(WTW)_1]a = e0b3 β )h(b), and can be solved for the coefficients. Model construction proceeds as indicated for Case 2 above. Hence, the model is unique only to within an arbitrary linear function. This analysis highlights a particular advantage of the current formulation, where the separate system (6.21) is derived for the two parameters d1 and d2. The specific mathe matical conditions for nonuniqueness in the calculated model can be determined simply by examining the equation det A = 0. Moreover, as indicated above, model construction can stifi proceed in these situations. In contrast, the alternative development of Parker et al. (1987) is very ambiguous regarding both of these issues. 6.3 Interpolation via the smoothest model The theory developed in the previous section is clearly illustrated by the problem of constructing a smooth interpolation of a set of discrete samples. This is a common problem in many branches of geophysics. The particular application considered here consists of calculating a continuous depth profile of a refracting interface from a set of point estimates of the refractor depth and dip. These point estimates can be derived by analyzing the first break traveltimes of a seismic experiment by various techniques. There are many such techniques available to the practicing interpreter, and these are not examined individually here. However, the methods typically employ the assumption that the critically refracting horizon is plane or nearly plane. Under these circumstances, it is logical to construct an interface depth proffle that possesses minimum curvature. Thus, the final model for the interface adheres closely to the prior theoretical assumptions used for inferring its depth. Many shallow seismic refraction projects are undertaken in conjunction with a driffing program (e.g., Hawkins, 1961; HasselstrΓΆm, 1969; Hatherly and Neville, 1986; Schwarz, 1990). Accurate depths to an interface are determined at the drillholes and the refraction 134 method is then used to extend this information between or beyond the wells. In some cases, outcrops of the target horizon provide additional geologic control on both depth and dip (Kilty et at, 1986). When constructing a proffle of the interface, the depth and dip information from all available sources should be integrated together. Extra data from dri]lholes, trenching, outcrops, etc. then contribute directly to the solution, rather than merely providing βtiepointsβ to the profile determined from refraction data. This goal can be readily achieved if the problem of caiculating the refractor proffle is posed as an interpolation issue. In this case, the data weighting matrix W assumes an important role. This matrix can be used to nondimensionaJize different data types, and to emphasize those particular data deemed more important in the inversion. In contrast, there is no provision in the GRM depth conversion technique for including additional geological or geophysical data with variable weighting. 6.3.1 General theory The complete set of N measured data is divided into two distinct subsets: n depth samples and N β n dip samples, where 0 m N. These data are observed at abscissae aj that are ordered as follows: Depth: a < X2 < < < x Dip: a Xn+1 < Xn+2 < < XN1 < XN b. The abscissae associated with the dip samples do not necessarily coincide with those of the depth samples. In particular, note that the extreme cases n = N (no dip samples) and n = 0 (no depth samples) are allowed. Expressions for the kernel functions pj(x; ci, c2) are now derived. The depth data are treated first. For j = 1,2,.. . , n, a predicted datum is considered to be a sample from a 135 function e(x): eTd e(xj) = The th component of the kernel function vector g(x) is a Dirac delta function centered at the sampling point 2j : gj(z) = S(x β z,). Hence, the singly and doubly integrated kernels are h(x) = k1(x) An expression for the more general kernel function pj(z; c, c2) is then obtained from equation (6.10): p3(x;cl,c2) = R(xβx)βR(x βC1)+(Xj βci)H(x-βc2). (6.24) The clip data are examined next. A simplification arises if these data are redefined as slope, rather than dip samples. The transformation is easily effected by calculating the tangent of each dip angle. From the fundamental theorem of calculus, the difference in the slope of the function e(x) between two locations z1 and C2 is given by Zj & eβ(xj) β eβ(c2) = Β£2 eβ(z) dx = j [H(z β c2) β H(x β xa)] eβ(x) dx. Thus, for i = n + 1, n + 2, .. . , N, a predicted slope datum is a functional of the second derivative of e(x): eβ(x) = d2 + jbpi( z;cl,c2)eII(x) dx, where d2 = eβ(c2) and the kernel function is defined by pj(x;cl,c2) = H(xβc)βH(x-βxj). (6.25) 136 Comparison of the above two expressions with equations (6.9) and (6.10) indicates that h(b) = 0 and k(b) = β1 for the index range j = n + 1, m + 2,. . . , N. Hence, if the dataset consists entirely of slope samples (n = 0), then the interpolation problem is an example of the Case 1 phenomenon described in the previous section (i.e., the constructed interpolant e(x) is unique only to within an additive constant). If the weighting function w(x) equals unity, then explicit formulae for the elements of the inner product matrix F(ci, c2) can be derived by substituting the above kernel functions into (6.16) and integrating. Arbitrary nonuniform weighting requires that the integrals be evaluated by numerical techniques. The examples presented below adopt uniform weighting. For Cl c2, there are 15 distinct formulae for the matrix elements Iβjj, corresponding to various locations of the coordinates z and x relative to c1 and C2. These equations are omitted here for brevity. If Cl > c2, then a transformation of the independent variable via = a + b β x and a reindexing of the data allow the previously established formulae to be used. As indicated previously, it is possible to solve for a, d1, and d2 simultaneously. This is the preferred method of solution for situations with a large amount of data. It eliminates the need to invert F(t) in order to obtain the elements of the matrix A and the vector b as in (6.22). However, for small scale problems like the following examples, calculation of F(t) is not computationally burdensome. Singular value decomposition (SVD) is used here to compute the inverse. An advantageous feature of this technique is that it allows the singular vectors associated with very small eigenvalues to be winnowed from the solution simply by truncating the SVD sum. Small eigenvalues arise when the matrix F(Au) is nearly singular, and contribute to a large variance in the computed solution. Let P(ji) denote the th element of the matrix {F(ci,c2)+1(WTW)_1] . Then, define three sums as follows: n ii N = > lββ(p), = x, = 3=1 j=rn+l 137 These quantities are all defined for the full index range i = 1,2,. . . , N. After some algebraic manipulation, the elements of the matrix A and the vector b in equation (6.21) reduce to the simple forms all P, a12 = (xj β ci )P + i=l 2=1 i=n+l n n N a21 (Q+Ri), a22 = (xβc1)(Qi+R) + (Q+Rj), i=l i=l = b2 = Solution for d and then the coefficient vector a is now possible. The problem is com pleted by constructing the function e(z) via equation (6.23). Expressions for the doubly integrated kernel functions in (6.23) are obtained by integrating equations (6.24) and (6.25) for p,(z; Cl, c2). There are five distinct cases that must be examined, each corresponding to different location of the coordinate x relative to the two abscissae Cl and C2. 6.3.2 Numerical example Figures 6.1 through 6.4 present a synthetic example from shallow refraction seismology. The top panel in each figure depicts the elevation profile of a near surface alluvium/bedrock interface (dashed line). A small channel is located adjacent to a larger and broader anticine. Elevations are referenced to an arbitrary horizontal datum plane. Nonuniformly spaced elevation and slope samples from the interface are indicated by asterisks. Initially, these are considered to be error free. If accurate estimates of interface depth and dip are available (say, from a shallow bore hole) then these values can be used to construct a smoothest elevation model. Figure 6.1 depicts this situation. The prescribed elevation and slope values are indicated by open 138 β02 0 10 20 40 50 β 0.20 : 50 x(m) Fig. 6.1. A smoothest refractor model e(x) constructed from six error-free depth samples (indicated by asterisks). Dashed lines (- - -) refer to the true earth model and solid lines (β) refer to the constructed model. The independent values of the refractor elevation and slope are taken from the true earth model and are indicated by small squares at abscissa ci = C2 = 23.75 m. The 12 norm of the second derivative of the constructed model is 0.228 m_h/2. 6.0 0 10 02 40 50 139 0.2 0.0 β02 0 10 20 0 40 50 0.10 _v βo.io 0 10 20 30 40 50 x(m) Fig. 6.2. Same as Fig. 6.1 except that optimum values of the refractor elevation and slope are calculated at Cl = c2 = 23.75 m. The 12 norm of the second derivative is reduced to II eβ = 0.129 mβ2. 6.0 0 10 40 50 -I I I I 140 squares at horizontal coordinate ci = C2 = 23.75 m. The constructed model and its deriva tive (solid curves) pass through these points. However, the second derivative of the model possesses a jump discontinuity at this location. This discontinuity arises from the Heaviside step function with onset at a = c2 in equation (6.24). If the depth and dip values are spec- fled erroneously (i.e., unequal to the true earth model values) then the cusp in eβ(x) and the discontinuity in eβ(x) can be magnified (Aidridge et al., 1991). Figure 6.2 displays the smoothest elevation model constructed by using the optimum values of d1 and d2. There are obvious differences from the previous model on the flank of the anticine and on the flat portion between anticine and syncine. The second derivative of the elevation proffle is now a continuous, piecewise linear function. It is emphasized that the same model is obtained if the parameters d1 and d2 are optimized at any pair of horizontal positions ci and C2 along the proffle (including outside of the range encompassed by the six elevation samples). The 12 norm of the second derivative is reduced by optimizing the parameters, as expected. Aidridge et at (1991) discuss the relation of these interpolants to cubic splines. In particular, the constructed elevation profile in Figure 6.2 is a cubic spline on [a, bJ, whereas e(z) in Figure 6.1 consists of two cubic splines joined end-to-end at the coordinate x = C2. Inclusion of slope samples, as in Figure 6.3, results in a more accurate reconstruction of the refracting interface and its first derivative. However, continuity of eβ(x) is now sacrificed. As predicted by equation (6.25), a jump discontinuity is introduced into the second derivative at each abscissa x associated with a slope datum. Hence, the constructed elevation proffle is no longer a cubic spline on [a, bj. In actual practice, the calculated interface depth and dip samples will contain some random error. This situation is simulated in Figure 6.4. Random numbers drawn from a uniform probability distribution on Β±0.50 m (standard deviation = 0.29 m) are added to the accurate elevations. Also, each slope sample is perturbed by an amount corresponding to a uniformly distributed random dip angle on Β±5Β° (standard deviation = 2.9Β°). A model that is an exact fit to these error contaminated samples would contain spurious structure induced by the noise. Hence, a model is constructed that reproduces the erroneous data 141 0.20 β. 0.10 V β0.10 50 Fig. 6.3. A smoothest refractor model constructed from six depth samples and six dip samples. All data are error-free. Optimum values of the refractor elevation and slope are calculated at Cl = C2 = 23.75 m. eβ 0.354 mβ2. 6.0 βp 0 10 20 β0.2 40 0 10 20 40 50 0 10 20 30 40 x(m) 142 7.0 16.0 50 0.2 nN -β-β 00 4, β02 -4 . 0.10 N - __ __ v β0.10 0 β10 20 30 40 50 x(m) Fig. 6.4. A smoothest refractor model constructed from error contaminated data. Depth and dip samples are located at the same abscissae as in Fig. 6.3. Optimum values of the refractor elevation and slope are calculated at c = C2 = 23.75 m. J= 0.259 mβ2. 0 tO 20 50 40 0 tO 20 SO 40 50 β 1β - β’1 β 143 approximately, rather than exactly. The tradeoff parameter ,u and the weighting matrix W can be used advantageously for this purpose. Since a constructed smoothest model depends on the value of the tradeoff parameter, predicted data generated by this model will also be a function of p. The notation ez,,.d(p) is used to explicitly denote this dependence. If the weighting matrix is diagonal with elements W:, i = 1,2,. . . , N, then the degree of misfit between observed and predicted data is quantified by Edepth(JL) j,2 {e β erd(p)]2, 6slope(p) N β W [ebs β e()j2. - jfl+1 A value for the tradeoff parameter is sought such that these misfits are approximated by depth(P) Eslope(P) Nβn Note that, for uniform weighting, the right hand side of each expression above is proportional to the rms value of the standard deviations of the data errors. Selection of a search procedure for p is largely a matter of personal preference. For small scale problems like the present example, a few trial runs with guessed values of p allow a sufficiently precise estimate to be made. Alternately, it is possible to implement an automated iterative scheme that converges upon the desired tradeoff parameter. Figure 6.4 displays a refractor elevation model constructed with the misfit values Edepth = 0.30 m and e3j0 = 0.0063. The slope data are weighted by 0.12 relative to the depth data. Hence, the misfit to the slope samples corresponds to one standard deviation of a uniform distribution of dip angles on Β±5Β°. Also, optimum values for the two parameters d1 and d2 are used. Interest in detecting the shallow channel may stem from various exploration 144 objectives: placer ore localization, groundwater accumulation, contaminant migration, etc. However, the figure suggests that the channel is just at the limit of detectability for the given horizontal sampling and error statistics. 6.4 Conclusion The conventional method for constructing the smoothest model is generalized in two important ways: 1) The additional information required to calculate the smoothest model may be specified anywhere within the interval of definition [a, bj. This improvement can be particularly useful if the best known value of the model (or its slope) is not located at an endpoint. 2) Madmum smoothness can be obtained by calculating optimum values of the model and its slope directly frorri the observcd. data. This technique is advisable if independent values of these parameters are unknown or are difficult to estimate accurately. It is evident that these techniques can be extended to model construction situations where the 12 norm of a higher derivative of the model is minimized. The modified data equation satisfied by the m derivative of the model m()(z) will entail ii parameters cij de fined at the abscissae cj as follows: dj m(1)(cj), i = 1,2,... ,n. Extremizing the relevant objective function leads to an n x n system of equations Ad = b for the unknown d. A unique solution exists when det A 0, and the model constructed with these parameters will possess a continuous th derivative (m(a) E Cβs) provided that the original kernel func tion vector g(x) E C (this assumes that the weighting function w(z) is also continuous). If the determinant vanishes, then the model that yields an absolute minimum of m(β) 12, subject to constraints provided by the data, can still be constructed but is only unique to within an additive polynomial of at most degree n β 1. Only the smallest model (n = 0) is unique in all situations. Application of this model construction formalism to the problem of calculating an in terface depth profile yields a useful tool for seismic refraction exploration. The smoothest model is the natural model to adopt in this situation because refraction traveltime inversion 145 methods assume, either explicitly or implicitly, that local interface curvature is negligable. The resulting depth profile and its derivatives are mathematically well defined throughout the horizontal range (a, b) and thus can be used for other forward modeling purposes. The method is also flexible, as evidenced by the following specific advantages compared with depth determination via the GRM: 1) Additional depth and clip data arising from a variety of geological, geophysical, or engi neering techniques are readily incorporated into the computation. 2) Variable weighting of all data is easily achieved. 3) An adjustable misfit to error contaminated data is possible. Numerical examples presented here and in Aldridge et al. (1991) demonstrate the feasi bility of the technique for refractor depth profile construction. Obviously, there is stifi much to learn about issues such as weighting of different data types, resolvability of small scale features, optimum horizontal sampling, etc. However, the algorithm presented here provides a flexible tool for investigating these and related phenomena. 146 CHAPTER 7 REFRACTOR IMAGING USING AN AUTOMATED WAVEFRONT RECONSTRUCTION METHOD 7.1 Introduction The Wavefront Method is one of the earliest of the many techniques for interpreting re fraction arrival times. In 1930, Thornburgh demonstrated that subsurface wavefronts could be reconstructed from surface arrival times by applying Huygensβ principle in reverse. Sub sequently, Hagedoorn (1959) elucidated an imaging condition for delineating a refracting horizon. First, two oppositely propagating wavefront systems are reconstructed from the arrival times recorded on a forward and reverse spread, respectively. Then, pairs of these subsurface wavefronts intersect on or slightly below the refracting interface when the sum of their times equal.s the known reciprocal time (the shot-to-shot traveltime). This imag ing principle yields the correct spatial locus of a critically refracting horizon if the earth consists of constant velocity layers βbounded by plane dipping interfaces. However, several investigators have demonstrated that the imaging condition is reasonably accurate even if the measured arrival times are due to diving rays, rather than true critically refracted rays (Hagedoorn, 1959; Rockwell, 1967; Schenck, 1967; Hill, 1987). Diving rays may arise from nonpiane structure on the refracting interface, or a velocity gradient within the underlying medium. Extensive application of the wavefront method has been limited by two factors: i) la borious graphical techniques are required to construct the subsurface wavefront loci, and ii) detailed knowledge of the near surface velocity structure is necessary. This study directly ad dresses the first of these two issues. Instead of defining the wavefronts by a tedious graphical application of Huygensβ principle (e.g., Rockwell, 1967), a finite-difference computer algo rithm is used to downward continue surface arrival times through a specified velocity field. The algorithm is rapid and accurate, and is capable of handling a heterogeneous velocity structure. 147 Recently, Hill (1987) downward continued refracted waveforms to obtain a two dimen sional image of shallow structure. Although the technique presented here works only with arrival times, the goal is identical. The advantage of this approach resides in its compu tational simplicity. Since the propagation algorithm operates directly in the space-time domain, no transformations of the recorded wavefield, with attendent concerns about sam pling adequacy (Clayton and McMechan, 1981; Hill, 1987) are necessary. Furthermore, true amplitude recording and processing of seismic traces are not required. However, prior pick ing of these traces to obtain the arrival times is necessary, and this may be a time consuming job in some situations. 7.2 Finite-difference traveltimes 7.2.1 Wavefront Construction Vidale (1988, 1990) has recently developed an algorithm for calculating the first arrival times of a seismic wave propagating through a two or three dimensional velocity structure. The velocity field is sampled on a uniformly spaced 2D or 3D grid; plane wave finite-difference operators are used to extrapolate the traveltimes from point to point throughout this grid. Calculations are initiated at a source point within the predefined velocity field. The al gorithm properly handles the various wave types that comprise first arrivals (body waves, head waves, and diffractions). Subsequent contouring of the computed traveltime field yields a visual impression of propagating wavefronts. The phrase wavefront construction is used here to refer to traveltime loci calculated in this manner. Figure 7.1 depicts the subsurface wavefront systems generated by a sequence of shots buried in an earth model with undulat ing surface and refractor topography. The direct wave through the overburden is the initial arrival near each shot location. Beyond the crossover distance, the wave refracted by the higher velocity bedrock arrives first. The traveltimes recorded along the nonplane surface of the model are accurately computed by assigning a P-wave velocity to the uppermost layer equal to the speed of sound in air (β-i 350 m/s). Wavefront contours are then suppressed in de pt h (m ) 20 10 de pt h (m ) 0 30 20 10 30 0 aβ 0 30 0 de pt h (m ) 20 10 0 30 0 de pt h (m ) de pt h (m ) 20 10 0 30 20 10 0 β C n C1 2 e 4 - . , β . b .β’ c - ,- .- C C D I I . C) β e 4 I) Γ· : z β’ Q CD β - CD )-j p, I- , o C β’ - D - 0 0 CD CD p, )- CD CD . < CI ) CD q CD Cβ ) C n _ I-+ -. CD , 0 C.) $ CD CD CD CD II I- β. CI :, Cl ) E o β - - 0 N 0 a β (I) I- . 0 La (7 β 0β 0 -1 0 I 0 0 I 0 0 - 2 01 149 0 t/0 0 horizontal position (m) Fig. 7.2. Surface arrival time curves for the five wavefront systems depicted in Fig. 7.1. this region for visual clarity. The surface arrival time curves displayed in Figure 7.2 illustrate that nonpiane topography has a complicating effect on an interpretation. Vidaleβs wavefront construction algorithm has been altered in two important ways in order to improve its suitability for the shallow refraction problem. First, the traveltime calculations are initiated from a spatially extended source, rather than a point source. In the 2D case, source activation times are specified on and inside a rectangular region located within the velocity field; arrival times at grid points outside this rectangle are generated by the normal working of the algorithm. Since wavefronts are strongly curved in the immediate vicinity, of a point source, use of the plane wave finite-difference operators will yield inaccu rate traveltimes in this region. Moreover, these inaccuracies will be propagated to greater distances, where the plane wave extrapolators are locally valid. In order to avoid this prob lem, near source traveltimes are calculated via mathematically exact formulae appropriate 0 25 50 75 100 150 for either a constant or linear velocity field. Although more complicated velocity distribu tions can be considered, these particular velocity functions provide sufficient flexibility for many traveltime computation problems. Second, the mathematical form of the traveltime extrapolation operator is modified in those cases where there is a large velocity increase across a grid cell. This situation is rela tively common in the shallow refraction environment. The interface between unconsolidated overburden and consolidated bedrock, or between saturated and unsaturated alluvium, often represents a sharp velocity increase. In these cases, as the following analysis indicates, the conventional traveltime extrapolation formula may fail. Figure 7.3a depicts a system of plane wavefronts propagating across a square grid cell with side length h. The arrival time at the corner numbered 4 must be calculated from the known arrival times t1, t2, and t3 at the other three corners of the cell. Assuming a plane wave advancing with a constant slowness s, this time is given by t4 β ti + (/h cos where the angle S describes the ray direction relative to the cell diagonal. Simple geometric analysis yields I t3_t2 COSS = 1 β ___ _ Hence t4 = t1 + i/2(hs)2 β (t3 βt2) . (7.1) This expression is identical to Vidaleβs (1988) equation (3), which was derived by approx imating the partial derivatives in the 2D eikonal equation by finite differences and then solving algebraically for t4. The present derivation clearly reveals the underlying geometric assumption of plane wave propagation. If the argument of the square root in equation (7.1) becomes negative, then the plane wave extrapolation formula is obviously invalid. This may occur, for example, if there is a dramatic velocity increase across the cell (implying that the slowness s assigned to the cell 151 Fig. 7.3. Finite-difFerence traveltime extrapolation operators. (a) Locally plane wavefronts; 9 is the angle between the ray (heavy line) and the cell diagonal. (b) Circular wavefronts. In each case, the square grid cell has side length h and assigned slowness s. is quite small). In these cases, the arrival time t4 is calculated via the alternative formula = ruin {t1 + t2 + hs, t3 + hs}. (7.2) The geometric basis of equation (7.2) is ifiustrated in Figure 7.3b. In effect, the plane wavefront approximation is abandoned and Huygensβ principle is used directly to calculate the next traveltime. Although this computed time is not exactly correct, extensive numerical testing indicates that equation (7.2) is superior to the fix advocated by Vidale (1990) (i.e., if the argument of the square root becomes negative, take t4 = t1). a) ti t2 b) ti t2 h t3 t4t3 t4 \ \ 152 7.2.2 Wavefront reconstruction Figure 7.4 displays the forward and reverse wavefront systems generated by shooting over a shallow syncine. These are the wavefronts that give rise to the first arrival times observed on the surface. The finite-difference traveltime algorithm can now be used to recreate subsurface wavefronts from knowledge of the arrival times recorded at the surface. The source rectangle is placed at zero depth and is greatly elongated in the horizontal dimension. This line source is then activated sequentially (rather than simultaneously) with an initiation function T(z) derived from the recorded refraction arrival times T(z): Ts(X) = TR β T(X), (7.3) where TR is the reciprocal time. At source-receiver offsets X less than the crossover distance, phantom arrival times T(X) can be constructed from parallel traveltime curves recorded from distant shotpoints (Rockwell, 1967; Ackermaun et al., 1986). The line source generates a set of wavefronts radiating downward into the specified velocity field (Figure 7.5). The downward continuation velocity function v(z, z) is selected as a good approximation to the actual near surface velocity structure. Hence, within the overburden, the calculated wavefronts coincide with the emerging refracted wavefronts of Figure 7.4. Since the position of the refracting interface is initially unknown, the wavefronts are continued to greater depth using the known velocity field v(c, z). Rockwell (1967) referred to these traveltime loci as a βdirected wavefront systemβ; in this study, the phrase wavefront reconstruction is used to describe the process of creating an emergent wavefront system from recorded surface arrival times. 7.3 Refractor imaging Let tf(x, z) and tr(z, z) refer to the subsurface traveltime fields reconstructed from the forward and reverse arrival times, respectively. Then, according to Hagedoornβs (1959) imaging principle, the refracting interface is implicitly defined by the relation 153 2.5 I ci) 2.5 Fig. 7.4. Forward and reverse first arrival wavefronts for a shallow syndine model. Over burden velocity v = 1500 m/s, bedrock velocity v2 = 2500 rn/s. Grid cell size is 5 m and contour interval for wavefronts is 30 ms. 0 .5 0 1 1.5 2 .5 1 1.5 horizontal position (km) 2 154 I a) 2.5 I a) U) β’0 Fig. 7.5. Emergent wavefronts reconstructed from the surface arrival times recorded over the shallow syndine. Downward continuation velocity v(x, z) = 1500 rn/s. Grid cell size is 5 m and wavefront contour interval is 30 ms. 0 .5 0 1 1.5 2 .5 1 1.5 2 2.5 horizontal position (km) 155 tf(X,Z) + tr(x,z) TR. (7.4) Figure 7.6 graphically illustrates the superposition of the two reconstructed wavefront sys tems shown in the prior figure. This 2D array of superposed traveltimes is systematically searched to locate grid points where the imaging condition (7.4) is satisfied. If equation (7.4) does not hold on a grid point, linear interpolation between adajacent points is used to find the proper depth. The resulting depth locus z() (dashed line in Figure 7.6) is an accurate spatial image of the original refracting horizon, except near the edges of the input velocity field where the subsurface wavefronts are not reconstructed correctly. Figure 7.7 indicates that the technique is also capable of imaging anticinal structure The apex of the anticine is imaged slightly too deep because the refracted rays penetrate beneatk this stucture, rather tKan propagating along the undulating interface (e.g., Rage- doom, 1959, figures 2 and 3). Note that a. similar problem does not occur with the syndine, because there is a tendency for the diffracted ray to follow the interface in the presence of syncinal structure. The calculated locus for the refracting horizon depends on the reciprocal time TR and velocity field v(x, z) used for downward continuation of the surface arrival times. Variations in these quantitites from their correct values will induce variations in the depth and position of the refractor. It is relatively easy to assess the dependence of the refractor image on the value of the reciprocal time. The forward and reverse subsurface traveltime fields are added together and the result is contoured for various candidate βimaging timesβ. Figure 7.8a illustrates this situation for the buried syndline. If the imaging time used is less than or greater than the true reciprocal time, then the interface image is too shallow or too deep, respectively. This particular dependence upon the reciprocal time is the converse of that predicted by the classical wavefront method (Rockwell, 1967, p. 378). The difference arises from the method of reconstructing the subsurface wavefronts. Currently, finite-difference traveltime computations are initiated with the source function (7.3) and the algorithm is run forward 156 0 I U) 0 2.5 0 I 1) 0 .5 1 1.5 2 2.5 horizontal position (km) Fig. 7.6. (a) Superposition of the two reconstructed wavefront systems of Fig. 7.5. Dashed line is the locus satisfying the refractor imaging condition. (b) Comparsion of the true (solid) and imaged (dashed) refracting interfaces. 0 .5 1 1.5 2 I I I I b) I I 157 0 I 2.5 0 I .4) U, Cβ2 0 .5 1 1.5 2 2.5 horizontal position (1cm) Fig. 7.7. (a) Superposition of the forward and reverse reconstructed wavefronts (contour interval = 30 ms) over a shaJiow anticine. Dashed line is the refractor image. (b) Comparison of the true (solid) and imaged (dashed) refracting interfaces. U) .5 1 1.5 2 I I I I b) I I I 158 0 I ___ ___ ___ ___ ___ ___ Lt) C2 .5 1 1.5 2 2.5 0 I r - Li) e 2.5 horizontal position (km) Fig. 7.8. Dependence of the syncline locus on imaging time (top) and downward continuation velocity (bottom). Different images correspond to increments of 10 ms in imaging time and 100 rn/s in velocity, respectively. Heavy lines are the images corresponding to the correct values of TR and v. I I I __ _ _ _ _ __ _ TRβ6O ms 1TR+6O ms 0 .5 1 1.5 2 159 in time. Hence, subsurface wavefronts are labeled with times later than the surface source values. In contrast, the classical wavefront reconstruction methods label the subsurface wavefronts with times earlier than the surface measured times. In either case, the true position of the interface corresponds to an imaging time equal to TR. Quantifying the dependence of the refractor image on the downward continuation ve locity is more complicated. Forward and reverse subsurface wavefront systems must be reconstructed for each velocity function used in the analysis. Figure 7.8b displays a set of images of the shallow syncine calculated for various values of a constant downward continu ation velocity. If the velocity is less than or greater than the actual overburden velocity, then the interface image is too shallow or too deep, respectively. Moreover, a grossly incorrect continuation velocity distorts the shape of the interface structure. Hence, in common with many other seismic refraction interpretation techniques, accurate time-to-depth conversion with the wavefront method requires good knowledge of the overburden velocity distribution. This information can be obtained from uphole times, direct and reflected arrivals, shallow refractions, and borehole data. Finally, the accuracy of the solution depends on the reliability of the picked first arrival times. The ability of the method to resolve small scale features on the refracting horizon is also limited by the field geophone interval. These phenomena are analyzed by performing the inversion with noisy traveltime data sampled at an assumed geophone interval. Figure 7.9 displays the refractor image obtained by downward continuing error contaminated arrival times sampled every 25 m. Spatially correlated, normally distributed time errors (standard deviation 5 msec; correlation distance = 100 m) are added to the theoretically exact refraction picks. Appendix F describes the method used for computing correlated random noise. A cubic spline is then loosely fitted to the noisy arrival times and is used in equation (7.3) for the source initiation function. The experiment has detected the presence of the syndine, and its lateral position and depth are approximately correct. Long wavelength undulations on the refractor are artifacts of the spatially correlated noise, but are not unduly harmful to the structural interpretation. 160 0 I 4-β LI) 0 2.5 Fig. 7.9. Reconstructed wavefonts and syndine image formed from noisy traveltime data sampled at a geophone interval of 25 m. Grid cell size is 5 m and wavefront contour interval is 30 ms. Long wavelength undulations on the refractor arise from spatiafly correlated traveltime errors. 7.4 Refractor velocity estimation A particular advantage of the wavefront method is that the interface depth calculation is independent of the refractor velocity. Rather, the velocity of the substratum can be es timated after the position of the refracting horizon is determined. The distance between two points on the interface divided by the difference in the reconstructed wavefront times at these points is an estimate of the refractor velocity. This value is assigned to the midpoint of the two points for plotting purposes. In effect, the directional derivative of the subsurface traveltime field along the interface locus is computed by a centered finite-difference formula; the reciprocal of this value corresponds to the local velocity of the refractor. Either the for ward or reverse wavefront systems may be used for the computation. Figure 7.10 illustrates that the refractor velocity estimated in this manner possesses systematic errors related to .5 1 1.5 2 horizontal position (km) 161 CO (I) Fig. 7.10. Refractor velocity estimates for the syncine (top) and anticine (bottom) models. The spatial differencing interval used for the calculation is 200 m. Small amplitude, short wavelength oscillations are artifacts of the grid interval. Co cβi .75 1 1.25 1.5 1.75 .75 1 1.25 1.5 1.75 horizontal position (km) 162 the interface structure. However, for the two synthetic examples examined here, the inferred velocity values are everywhere within 3% (or Β± 75 m/s) of the correct values. 7.5 Field data example The interface imaging procedure is tested with a shallow refraction dataset acquired at the archeological site of Phalasarna in western Crete. Hadjidaki (1988) discusses the historical and archeological significance of this site and also gives a detailed description of the surface and near subsurface conditions. Forward and reverse refraction proffles were recorded along an inline spread of 18 geophones (geophone interval = 0.5 m, near source offset 0.5 m) during the summer of 1989. The data acquisition system consisted of a portable signal stacking seismograph with a hammer energy source. First arrival time picks were made on stacked traces in order to reduce random errors induced by ambient noise. The arrival times observed at reciprocal source positions at opposite ends of the spread differ slightly (16.40 msec vs. 16.70 msec on the forward and reverse profiles, respectively).. Since a. successful inversion requires consistency in the measured reciprocal times, reciprocal time corrections (Hatherly, 1982) are applied to the picked arrival times. A constant time shift is added to the raw time picks on each source gather in order to adjust the observed reciprocal times to the average value of 16.55 msec. Figure 7.11 displays the 18 first break picks recorded on the forward and reverse proffles after application of these reciprocal time corrections. A preliminary interpretation of the plotted traveltime curves identifies the direct and refracted branches. Overburden velocity, determined from the slopes of the direct arrival segments, exhibits a weak lateral variation (βββ 8%) over the 9 m spread length. This information is used to construct a near surface velocity function for subsequent downward continuation of the refracted arrival times. A cubic spline is fitted to the 16 refraction picks on each spread and is extrapolated to zero offset as a straight line. These curves are then used in equation (7.3) to calculate the source initiation functions required for the wavefront reconstruction algorithm. 163 horizontal position (m) Fig. 7.11. Shallow refraction traveltime data acquired at the dry harbor of Phalasarna, Crete. First break picks are indicated by triangles, and interpreted arrival time branches by smooth curves. Figure 7.12a depicts the subsurface wavefronts generated by downward continuing the refracted arrival times through a near surface velocity field given by v(x) = 252 β 2.11 (x in m and v in m/s). A shallow undulating interface is then imaged using the corrected reciprocal time TR = 16.55 msec. The refractor velocity estimate (Figure 7.12b) exhibits two distinct zones: i) abrupt variations about 800 rn/s on the left, and ii) a low velocity zone slower than 800 rn/s on the right. The validity of various refractor velocity functions can be tested by using the wavefront construction algorithm to compare predicted traveltimes with the observed traveltimes. A uniform refractor velocity of 800 m/s leads to unacceptably large differences (Figure 7.13b). When a low velocity zone is introduced into the refractor, good agreement is obtained (Figure 7.13c). The effect of the low velocity zone is evident in the plotted forward and reverse wavefronts of Figure 7.14 between 6 m and 8 m horizontal position. Note that 0 0 1 2 3 4 5 6 7 8 9 r - . - . Ct, CD I- . 0 β-4 β - i. - - CD . CD CD Γ§β a q U ) CD 0 . . I 4 .9 CD - - β % β , o CD a. ,. 0 o j c 4 β’ 4 -β O 0 β β 0 1 -s U ) 0 I- . o C D 0 U ) - I . v e lo ci ty (m /s ) 20 0 40 0 80 0 80 0 10 00 0 2 de pt h (m ) 1 0 0 CA ) Ci ) 0 β . N 0 0 (1) c- I 0 Q1 Q1 - 2 - 2 c o Cl) -4-) C) 0 G.) CD β-4 a) 0 CD β-I Vi β-4 F 0 0 0 β-I 0 0 CD 0 0 CD 0 0 0 0 CQ0 165 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 horizontal position (m) Fig. 7.13. Comparison of observed and predicted first break times for two refractor velocity functions. (a) Constant velocity (dashed) equals 800 rn/s and variable velocity (solid) in cludes a low velocity zone between 6 and 8 m. (b) Predicted traveltimes calculated with the constant refractor velocity. (c) Predicted traveltimes calculated with the variable refractor velocity. .4-, 166 Fig. 7.14. Subsurface wavefronts (contour interval = 0.5 ms) constructed from the earth model with the laterally varying refractor velocity. Grid cell size is 0.02 m. Note the effect of the low velocity zone between 6 and 8 m. 0 .4-, C) 0 1 0 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 horizontal position (m) 167 the strong variations in refractor velocity displayed on the left in Figure 7.12b need not be incorporated into the model to obtain an adequate fit to the measured arrival times. Further adjustment of the refractor velocity to obtain a closer fit is probably unwarranted. Although the recovered model of Figure 7.14 generates acceptable predicted traveltimes, it cannot be stated with certainty that this is the correct earth model. Other interpretations of the observed traveltime data, incorporating multiple layers or lateral changes in structure and/or velocity, are possible. Since the refraction dataset does not include arrival times recorded from far offset shotpoints, it is not possible to distinguish between these alternatives (Ackermann et al., 1986). However, the inferred model is consistent with known subsurface information from the vicinity of the refraction proffle. Archeological treadling conducted about 17 m distant encountered dipping sandstone bedrock at approximately 1.8 m depth (Hadjidaki, 1988). Various earthen and gravel layers overlie the bedrock. A porous, aerated sandstone might have a P-wave velocity as low as β.β 800 rn/s. Hence, the preliminary interpretation is that the interface imaged in Figure 7.14 is the upper surface of the sandstone bedrock. The very low velocity between 6 m and 8 m may be a zone of more extensive weathering, fracturing, or aeration. Alternately, it is possible that one of the overlying shallow gravel layers has been imaged. 7.6 Conclusion The essential requirements for reconstructing shallow refracted wavefronts are: 1) arrival times T(X) from a given marker horizon recorded β’ (or phantomed) on a forward and reverse spread, 2) a reciprocal time TR, 3) a near surface velocity function v(x, z). A simple modification of Vidaleβs finite-difference traveltime algorithm then allows the rapid calculation of the subsurface wavefront systems that give rise to the recorded arrival times. 168 Although the synthetic examples described here use a uniform near surface velocity, down ward continuation through a varying velocity field is also possible with no increase in com putation time. The buried refracting horizon is delineated in a subsequent step by applying Hagedoornβs imaging principle. No prior assumption regarding the refractor velocity is re quired. Rather, the velocity of the substratum can be estimated by calculating the directional derivative of the reconstructed wavefront systems along the imaged interface. Picking of first arrival times and assignment of these picks to specific refractors are necessary in this method. The final locus for the refracting interface is sensitive to errors - in the picked times, as well as to an incorrect choice of the reciprocal time and velocity field. However, since the technique is not computationafly intensive, it is possible to assess the magnitude of the position and depth uncertainty by performing the inversion repeti tively. The forward modeling capabffities of the finite-difference traveltime algorithm can also be used to quickly generate predicted arrival times from the inferred subsurface model. Comparison of these times with the observed data is a powerful method of establishing the significance of various features of the recovered earth model. Finally, two specific problem areas with the automated wavefront reconstruction method have been identified that merit further research: i) downward continuation of traveltimes recorded along a nonplane surface, and ii) correction of the refractor velocity function for the effects of structurally induced errors. Although a fully automated solution to these problems is not yet available, this should not prevent the immediate application of the method in shallow seismic refraction exploration. 169 CHAPTER 8 TWO DIMENSIONAL TOMOGRAPHIC INVERSION WITH FINITE-DIFFERENCE TRAVELTIMES 8.1 Introduction Curved ray traveltime tomography was originally developed by Bois et al. (1972) for the purpose of estimating the seismic velocity distribution between two boreholes. Following their seminal work, several investigators advanced the technology of tomographic imaging with curved raypaths (Lytle and Dines, 1980; Bishop et aL, 1985; McMechan et aL, 1987; Bregman et al., 1989; White, 1989). Curved ray methods are necessary for accurately reconstructing the velocity field in highly refractive media. Although straight ray techniques are adequate in media with relatively small velocity variations, it is difficult to decide when the simplifying assumption of straight line raypaths becomes invalid. Hence, there is a compeffing reason to use curved ray methods in all situations: they are based on a more accurate model of wave propagation through variable velocity media. Traveltime tomography is a nonlinear inverse problem that can be solved by local lin earization and iteration. Since the velocity field is updated on each cycle of the tomographic imaging procedure, rays have to be retraced between all source-receiver pairs. This raytrac ing constitutes a large part of the computational cost of curved ray tomographic inversion. Two-point raytracing between a specific source and receiver is an iterative process, and may encounter difficulties due to shadow zones and multipathing. In this study, the problems associated with conventional raytracing are circumvented by calculating traveltimes to all points of a two dimensional slowness field with a rapid finite-difference algorithm (Vidale, 1988). Raypaths are then generated by following the steepest descent direction through the computed traveltimes from each receiver back to the source. This method yields the raypaths of all wave types that comprise first arrivals (body waves, head waves, and diffrac tions). Moreover, since arrival times are calculated throughout the slowness field, arbitrary recording geometries are easily accommodated. 170 In addition to rapid and accurate forward modeling, tomographic inversion requires the solution of a system of linear algebraic equations to obtain the improved velocity field. This solution may exhibit erratic and unphysical behavior due to noise in the observed traveltime data and/or ill-conditioning in the equations. Hence, some form of regularization is often used to stabilize the inversion. For example, Lytle and Dines (1980) introduce Laplacian smoothing into the system when calculating a perturbation to a slowness model. Bishop et al. (1985), Bregman et al. (1989), and White (1989) limit the size of the model perturba tion by using the damped least squares method. Macrides et al. (1988) impose inequality constraints on a perturbation calculated via an ART algorithm. The present approach is to apply linear equality constraints directly to the slowness model, rather than to a model perturbation, on each iteration of the inversion procedure. In addition to improving the mathematical conditioning of the system, the constraint equations allow the introduction of a priori geological or geophysical knowledge about the model into the iuversioa. In particu lar, the constraints may arise from a desire to impose a preferred character, like flatness or smoothness, on the slowness solution. Alternately, one may seek a model that is close, in some quantitative sense, to a prescribed base model. Inclusion of these constraint equations restricts the nonuniqueness that is common in realistic tomographic inverse problems. The nonlinear tomographic inversion procedure described in this chapter consists of four basic steps: 1) calculation of first arrival traveltimes from each source location to all points of a gridded slowness field, 2) generation of raypaths between all source-receiver pairs, 3) solution of a large and sparse system of linear equations for a perturbation to the existing slowness model, 4) updating and (optionally) smoothing the slowness model. This four-step process is initiated with an estimate of the true slowness function, and is repeated until an acceptable match is obtained between observed and calculated traveltimes. 171 The initial estimate is usually a uniform slowness field. Subsequent sections describe these steps in more detail, and demonstrate the inversion procedure using synthetic traveltime data from simulated VSP and crosswell experiments. 8.2 General theory and model representation The traveltime of a seismic wave propagating through a slowness field s is given by the path integral t(s) = f sdl, (8.1)r(s) where r(s) is the raypath connecting source and receiver, and dl is an incremental path length. Since the raypath locus depends on the slowness, the traveltiine t(s) is a nonlinear functional of s. The problem can be linearized by considering the traveltime difference At t(s + As) β t(s), where As is a perturbation to the slowness field s. For a sufficiently small perturbation, the raypath F(s + As) is approximated by the original raypath F(s). Fermatβs Principle then implies that the traveltime through the perturbed slowness field can be evaluated by integration along the unperturbed raypath. The traveltime difference becomes At = J As dl. (8.2)T(s) Hence, a smafl difference in traveltime is linearly related to a small difference in slowness. The two dimensional slowness function s(x, z) is represented by a set of K square cells, each with a uniform slowness value mk (k = 1,2, . .. , K). Thus, within a cell, a raypath is a straight line segment. For a collection of raypaths, equation (8.2) is expressed as the matrix/vector product At = A(m) Am, (8.3) where A(m) is a matrix of raypath length segments within the square cells of the slowness model m. In the tomographic inversion problem, a model perturbation Am is sought such 172 that the improved model m + m approximates the true slowness model that generates the measured traveltime data. Hence, the traveltime difference vector it is given by t0b3 t,.d(m), (8.4) where t3 is a vector of observed arrival times, and t,.d(m) is a vector of predicted tray eltimes computed from the known slowness model m. In principle, equation (8.3) can be solved for the required model perturbation. In practice, solution difficulties arise because the -raypath matrix is commonly nonsquare, large, sparse, and rank deficient. Moreover, since the observed traveltimes contain random errors, the system (8.3) may be inconsistent. In this case, an exact solution does not exist and a minimum misfit solution is usually sought. Finally, if the initial model m is a poor approximation to the true slowness, several iterations of the model updating procedure may be necessary before the inagnitude of the traveltime residual vector zt becomes acceptably small. 8.3 Forward modeling The first arrival times of a seismic wave propagating through a two dimensional ve locity structure are computed by Vidaleβs (1988) finite-difference scheme. This algorithm uses plane wavefront traveltime operators to extrapolate arrival times from point to point throughout a uniformly spaced grid. The method is rapid and accurate, and can be ap plied to a heterogeneous medium with moderate to strong velocity variations. Podvin and Lecomte (1991) describe improvements to the local traveltime extrapolators that allow mod els with very strong velocity contrasts to be examined. The traveltimes of all wave types that comprise first arrivals (body waves, head waves, and diffractions) are calculated. Reflections and other later arrivals are not included; this represents a limitation of the technique as currently formulated. Vidaleβs method is based on a centered finite-difference solution of the eikonal equation on each square cell of a gridded slowness field. Thus, the associated discretization error is second order in the grid cell size. An input slowness function s(x, z) is sampled on a 173 uniformly spaced two dimensional grid. If the grid interval is h, then the sampled slowness values are given by sj = sQc, zj), where ij = Xmin + (i β 1)h and Z = Zmjn + (j 1)h (with i 1,2,.. .,I and j = 1,2,... ,J). Since there are IJgrid points, there areK = (Iβl)(J--l) square cells. The slowness assigned to a particular cell is the arithmetic mean of the slowness values at the four bounding grid points: 1 mk = (sj + .s+i,j + si,j+i + S+1,j+1), (8.5) -where k i + (I β 1)(j β 1). With this indexing scheme, the elements of the vector constitute a row-ordered sequence of the two dimensional array of cell slowness values. An individual cell is referenced either by the coordinate indices of its upper left corner (ii), or by its sequential index (k). Calculations are initiated at a designated source point (,z3) (not necessarily coinci dent with a grid point) within the slowness model. Since wavefronts are strongly curved in the immediate vicinity of a point source, plane wavefront traveltime extrapolators are inappropriate in this region. Furthermore, near source inaccuracies are propagated to all greater distances. In order to mitigate these effects, the traveltimes in a near source rect angle are calculated via mathematically exact formulae derived from certain simple velocity distributions. Hence, consider the linear velocity function v(z,z) = v3 + a(x βx3)+a2(zβ z3), (8.6) where v3 is the velocity at the source location, and a and a are the horizontal and vertical components of the velocity gradient vector. Numerical values for these constants are obtained by performing a least squares fit to the velocity samples vj = l/sjj surrounding the source. The gradient components can be expressed in terms of the magnitude a and direction angle Γ§S (relative to vertical) of the gradient: a,, = a sin qβ and a = a cos q, where a = ../a2,, + a and tan 4 = a,,/a. If the frame of reference is rotated through the angle , then the velocity 174 field (8.6) is transformed into a one dimensional function. The spatial coordinates in the rotated system are designated by primes and are given by xl = x cos 4) β z sin 4), x SIn 4β + z cos 4β. In the rotated frame, the velocity function becomes v(zβ) = v3 + a(zβ β z). The arrival time at an arbitrary near source location (zβ, zβ) can then be calculated by standard 1D techniques. The result is I 21( i lβ2j( I\2 I l β β a LkX Xs) IZ Z31 LIx ,Z β Siflit ii ,a y 4v[v3+a(z βz3)J Wavefronts associated with this time field are eccentric circles with centers that are displaced along a stiaight line through the source point, in the direction of the velocity gradient vector. In the limit a β+ 0, it can be shown via LβHopitalβs rule that equation (8.7) reduces to the proper expression for a uniform velocity field (i.e., circular concentric wavefronts). The finite-difference algorithm calculates a traveltime tjj at every grid point of the slowness field. If a receiver is not located on a grid node, then an interpolator is needed to estimate the arrival time at the actual receiver position. Simple bilinear interpolation provides adequate accuracy. Hence, if a receiver with coordinates (x,., z7) is located within cell ij, then define the dimensionless quantities p = (z? β x)/h and q = (z,. β z)/h. The interpolated traveltime is given by t(xr,zr) = (1 βp)(l β q)t +p(l β q)tj,1+ q(1 βp)t,11 +pqt+i,+i. (8.8) If the four arrival times bounding cell ij are due to local plane wave propagation, then equation (8.8) is an exact expression for the traveltime at the interior point (x,., zr). In other cases, the interpolator (8.8) has accuracy 0(h2) (Dahlquist and BjΓΆrk, 1974, p. 319) and thus is consistent with the level of accuracy associated with the forward modeling scheme. 175 8.4 Ray generation Raypaths are generated by following the steepest descent direction through a computed traveltime field from each receiver back to the corresponding source point. This strategy was originally suggested by Vidale (1988), and has recently been implemented by Podvin and Lecomte (1991). The horizontal and vertical components of the traveltime gradient vector within cell ij are approximated by the centered finite-difference formulae Ot (t+1,+ t+,+i) β (tq + (8.9a) 2h (t,1 + β (t +t1,) (8 9b 2h Assignment of a constant traveltime gradient to a cell is compatible with the assumption of locally plane wavefronts used in the forward modeling algorithm. The steepest descent direction is opposite to the gradient direction of the traveltime field. Hence, within cell ij (or k), the steepest descent direction is defined by the angie β1 at/az = tan + ir (modulo 2w). (8.10) 8k is measured clockwise from the positive horizontal axis. For a fixed source, all raypaths that cross cell ij have this same orientation angle. The lengths of the raypath segments within the cell range from zero to a maximum of h. A representative situation for cell ij is depicted in Figure 8.1. The raypath enters the cell at point A on its right boundary with coordinates (ma, za). Depending on the value of the steepest descent angle assigned to the cell, the ray may exit on any of the remaining three sides or one of the four corners. The logic that selects one of these seven possibilities 176 xi xiii x β I I zi zj+1 z Fig. 8.1. Raypath (heavy line) traced through square cell ij of a 2D slowness model. Ray enters cell at point A, follows the local steepest descent direction across the cell, and exits at point B. is given in the first column of Table 8.1, where al and a2 are positive acute angles defined by β1 fZj+1 β Za\ β1 fZa Zj = tan h a2 = tan h These two angles are illustrated in Figure 8.1. After an exit option is selected, the coordinates of the exit point B are easily determined (columns 2 and 3 in Table 8.1). These coordinates (xl,, zb) constitute the entry point coordinates for the next cell that the ray crosses (columns 4 and 5 in Table 8.1). In the particular case displayed, the raypath enters cell i,j + 1 on its top boundary. Hence, a different logical scheme is required to extend the raypath across this next cell. A total of eight logical tables are necessary to handle all of the possibilities. B 177 Angle Range Exit Coordinates Next Cell Indices o < ir/2 i + 1 j + 1 ir/2 <8k <r β Xa + (Zj+i β za) cot 8k Zj+l i =lrβal Xi zi+i iβi j+1 1<9k.<lr+a2 xi zahtallOk iβi j xi Zj iβi fβi r + <8k <3ii-/2 x + (Zj β za) cot 0k i j β 1 3ir/2 9j <2ir z1 i + 1 j β 1 Table 8.1. Ray tracing logic for an entry point on the right side of a square grid cell. Column 1 gives possible ranges for the steepest descent angle 9k assigned to the cell. Angles al and a are defined in the text and ifiustrated in Fig. 8.1. Columns 2 and 3 give the horizontal and vertical coordinates of the ray exit point, respectively. (xa, za) are the ray entry point coordinates. Columns 4 and 5 give the indices of the next cell that the ray enters. These correspond to a raypath entering a cell on the right, bottom, left, and top sides, and the upper right, lower right, lower left, and upper left corners. The length of the raypath segment within cell ij is easily calculated once the entry and exit coordinates are known. Although tomographic inversion requires only the value of the segment length, the actual coordinates (xa, Za) and (zb, zb) are also retained in order to create raypath plots for diagnostic purposes. Ray tracing is initiated at each receiver position. If a receiver is located on a grid node or grid line, then an average of the steepest descent angles from the neighboring cells determines which cell the raypath enters first. Tracing then begins using the logic outlined above. However, if a receiver is positioned within a cell, a separate logical scheme generates the initial path segment to the enclosing cell boundary. The technique is similiar to that described above and is given in Appendix G. Iterative generation of path segments continues until the raypath arrives at the boundary of a defined near-source zone. The size and shape 178 of this region vary slightly depending on whether the point source resides on a grid node, a grid line, or within a cell (see Appendix G). The final portion of the raypath is then taken to be a straight line from the boundary point directly to the source position (c3,z3), regardless of the local values of the steepest descent angle. This ray termination procedure is designed to overcome difficulties associated with nonuniformity of the traveltime gradient vector in close proximity to the source. Figure 8.2a is a contour plot of a velocity model bounded by two vertical boreholes. A shallow low velocity anomaly overlies a dipping, higher velocity zone. The first arrival wavefronts from a surface source located between the boreholes are illustrated in Figure 8.2b. The wavefronts are retarded by the low velocity zone and advance more rapidly through the high velocity zone. In addition, Figure 8.2b displays the raypaths traced through this time field from 18 downhole receivers back to the source. The raypaths are orthogonal to the wavefronts, as expected. 8.5 Inversion mathematics As indicated previously, system (8.3) is typically ill-conditioned and inconsistent. Hence, an undamped least squares solution (m = (ATA)_lAT.t) may yield a model update vec tor with relatively large and unrealistic cell-to-cell variations in slowness. The conditioning of these equations can be improved by incorporating model constraint information into the inversion. Hence, the system (8.3) is augmented with sets of linear equality constraints on the updated model m + zm. In addition to mathematically stabilizing the inversion, these equality constraints allow the convenient introduction of a priori geological or geophysical knowledge (or bias) into the problem. In general, a set of linear constraint equations applied to the improved model m + m is written as the matrix/vector multiplication B(m + z.m) = b, (8.11) 179 250 horizontal position (m) Fig. 8.2. (a) Velocity model (contour interval = 100 m/s) with a horizontal low velocity zone and a dipping high velocity zone. Maximum velocity = 2800 rn/s; minimum velocity = 1633 rn/s. (b) Wavefronts (contour interval = 25 ms) and raypaths generated by a surface source. The raypaths are traced from 18 downhole receivers back to the surface source. 0 0 0 u0 0 250 500 r ci) 0 0 tb 500 180 where the coefficient matrix B and the right hand side vector b are prescribed. The subscript m (n = 1,2,. .. , N) is used to refer to a particular set of constraints. A perturbation Am that simultaneously satisfies, in the least squares sense, the linearized data equations and the model constraint equations is sought. Hence, the relevant objective function is (Am) = A Am β At 112 + B(m + Am) β b 112, (8.12) -where the scalars p (0 < < +oo) are adjustable tradeoff parameters that control the relative importance of the various terms. Extremizing with respect to Am yields the linear algebraic equations [ATA + p, BBn] Am ATAt + β B m). (8.13) For nonzero ji, the coefficient matrix in this expression is usually nonsingular. The re quired model perturbation Am can be obtained by solving (8.13) using standard techniques of numerical linear algebra. However, the coefficient matrix is large (K x K, where K is the number of slowness cells) and may be dense even if the original raypath and constraint ma trices are sparse. Hence, it is advantageous to seek a solution method that avoids explicitly forming the square matrices ATA and It is straightforward to demonstrate that (8.13) are the normal equations associated with the least squares solution of the rectangular system A At β Bim) Am = . (8.14) vhi(bN- BNm) 181 Algorithm LSQR (Paige and Saunders, 1982) is used to solve equation (8.14) directly for m. LSQR is an iterative solution technique for large and sparse systems of linear equations that is closely related to the conjugate gradient method. It is designed to seek the minimum norm least squares solution of a set of equations. Numerical studies of LSQR applied to the tomographic inversion problem indicate that it is both rapid and accurate (Nolet, 1985; Scales, 1987). A simple FORTRAN version of the LSQR algorithm is given by Nolet (1987). The dimensions of the coefficient matrix in (8.14) are (Ndata +N03) x K, where Ndati, and N3 are the number of data and constraint equations, respectively. Since this matrix may be large and sparse, a significant reduction in storage space is achieved by storing only the nonzero elements in a one dimensional array. The full index scheme described by Scales (1987) is used here to store and address the matrix elements. With this storage method, the sparse matrix/vector multiplications required by the LSQR algorithm are particularly simple to implement. - Finally, the correction to the slowness value at grid point ij is determined by averaging the slowness perturbations calculated for the four surrounding cells: (tXm_i + Amkj+1 + Amk_l + Amk), (8.15) where k = i + (I β 1)(j β1). Grid points located on the edges (corners) of the slowness model are updated by adding the average of the perturbations associated with the neighboring two (one) cells. After all grid points are updated, forward modeling of traveltimes for the next iteration of the inversion can proceed. Several investigators apply a spatial filter to the slowness field between tomographic iterations in order to suppress short wavelength variations in the computed values (Dines and Lytle, 1979; Radcliffe et al., 1984; Gersztenkorn and Scales, 1988). These variations may arise from noise in the traveltime data and/or instabilities in the inversion. Smoothing then serves to condition the slowness field for the next forward modeling step. Hence, an optional 9-point square smoother with prescribed weights Wi is included here. The smoother 182 is applied to the gridded values. The smoothed slowness value assigned to grid point ij is given by Sj = Wi Siβ1,jβ1 + W2 Si,j_1 + W3 Si+1,j_I W4 Siβ1,j + W5 S:j + W6 Si1,j tV7 iβ1,j+1 + W8 St,31 + W9 sii,ji , (8.16) where > = 1. Grid points residing on the edges of the model are smoothed by concep tuaJly extending the grid by one point with the local values. In the examples to follow, the ifiter weights are W5 4/12 and all other wj = 1/12. 8.6 Model constraint equations There is a wide variety of linear equality constraints that can be employed in the tomo graphic inversion problem. The particular type of constraint used in this study is charac terized by a two dimensional βoperatorβ or βfilterβ with five specified constantsci, C2, C3, C4, and C5. Figure 8.3 depicts the application of this operator to interior cell ij of a slowness model. Using the row-ordered indexing scheme, the kth component of the ifitered image Bm is given by Cl mj + C2 mk+I_l + C3 mk_l + C4 mk_I+l + C5 mj. Edge cells are handled by conceptually extending the slowness model beyond the defined region with the local cell slowness values. For example, application of the 5-point operator to the upper left corner cell (i = 1, j = 1 corresponding to k = 1) of the model m is via the formula ci m2 + C2 ml + (c3 + C4 + c5)ml. One such constraint is applied to each cell of the slowness model. Development of the matrix representation B for this set of equations is merely a matter of proper row and column 183 zi zj+1 - z C3 C5 Ci C2 Fig. 8.3. Schematic representation of the 5-point constraint operator applied to cell ij of a slowness model. The slowness values in five neighboring cells are multiplied by constants Cl through C5 and the products are summed. indexing. In general, the components of the right hand side vector b may differ, implying spatial variation in the value of the applied model constraints. The 5-point operator is an extremely simple and flexible mechanism for introducing various types of constraint information into tomograplilc imaging. Some examples include: Case 1: C5 = 1 and all other c = 0. An individual cell slowness is not constrained by its immediate neighbors. Rather, inversion produces a model that is closest, in the least squares sense, to a given target model b. Case 2: c = 1/2h, C3 β1/2h and all other c1 = 0. This operator is a centered finite difference approximation to the horizontal derivative of the slowness field. A similar approximation to the vertical derivative is achieved by setting c2 = 1/2h, C4 = β1/2h xi xi+i x I I C4 184 and all other ci = 0. These operators introduce first-difference regularization, or flatten ing, into the inversion. Case 3: ci = 1/h2, C5 = β2/h2,c = 1/h2 and all other cj = 0. This operator, as well as its vertical counterpart, introduces second-difference regularization, or smoothing, into the inversion. Case 4: c = C2 = C3 = C4 = 1/h2 and C5 = β4/h2. This operator incorporates Laplacian smoothing into the inversion. These constraints are inherently local in character; they consist of linear relations between adjacent cell slowness values. Obviously, an operator with larger spatial extent (say, a 9- point square ifiter) would provide more options. However, the current 5-point operator offers a reasonable compromise among the competing issues of flexibility, accuracy, and ease of implementation. - Combinations of constraints applied by the simple 5-point operator can also be con sidered. For example, the next section illustrates situations where both horizontal and vertical first-difference regulaΓ±zation is applied (i.e., N 2). The first system of constraint equations corresponds to c = 1/2h and C3 = β1/2h, while the second set corresponds to C2 = 1/2h and C4 = β1/2h. In each case, the right hand side vector b is set equal to 0. 8.7 Synthetic examples The examples presented in this section demonstrate the ability of the tomographic in- version procedure to image a smoothly varying velocity field. Figure 8.2a displays the 500 m x 500 m velocity model used to generate synthetic traveltime data. Since the grid interval is h = 5 m, there are IJ = 10201 grid points used for the forward modeling and K = 10000 square cells used in the inversion. The data acquisition geometry used for the first example simulates a double-well VSP plus crosswell experiment. Nine surface sources, located between two vertical boreholes, are spaced 50 m apart. Traveltimes are recorded by 18 borehole receivers (9 per well; 50 m 185 separation) from each source position. In addition, the 9 downhole receivers in the right well record traveltimes from 9 sources symmetrically placed in the left well. The complete set of 243 raypaths linking all source-receiver pairs is illustrated in Figure 8.4. Note that these first arrival raypaths tend to avoid the low velocity zone, resulting in a region of reduced ray coverage. Furthermore, no raypaths penetrate below 450 m depth. Contoured velocity tomograms obtained by inverting the combined VSP and crosswell traveltimes are displayed in Figure 8.5. As indicated above, both horizontal and vertical first-difference constraints are imposed in the iterative inversions. Additionally, in Figure 8.5b, the slowness values of the cells adjacent to the two boreholes are constrained by the true slowness function. In an actual field experiment, this information may be available from borehole velocity logs. Each of these reconstructions produces an rms traveltirne error of β-i 0.5 ms, which is about 0.25% of the rms value of the synthetic traveltimes (200.7 rns). The main features of the true velocity model are recovered by both of the inversions depicted in Figure 8.5. The location and amplitude of the dipping high velocity anomaly are approximately correct. Also, the shallow low velocity zone has been detected and correctly positioned, although the actual velocity value at its center is about 90 rn/s too high. This effect is associated with the reduced raypath density in this area of the model. Interestingly, inclusion of the borehole velocity constraints does not yield a dramatic improvement. The principal difference between the two reconstructions appears in the region below 400 m depth, where raypath coverage is negligible or nonexistent. Six (Figure 8.5a) and nine (Figure 8.5b) iterations are required to reduce the initial rms traveltime misfit to βββ 0.5 ms, with most of the improvement actually occurring on the first iteration. As stated above, the conjugate gradient solver LSQR is also an iterative algorithm. Theoretically, LSQR requires at most K iterations to converge to the solution of a system with K unknowns (assuming exact arithmetic can be performed). For the tomographic inversions described in this section, an acceptable solution is obtained with about 10 iterations in LSQR, which is three orders of magnitude less than the theoretical value K = 10000. This results in an appreciable saving in computation time. 186 Fig. 8.4. 243 raypaths traced through the velocity field of Fig. 8.2a. These rays link all source-receiver pairs of a combined double-well VSP and crosswell experiment. 0 β’4)O2 a, C 0 LOO 250 horizontal position (m) 500 187 500 250 500 horizontal position (m) Fig. 8.5. Reconstructed velocity models (contour interval = 100 m/s) obtained by inverting the combined VSP and crosswell traveltirnes. (a) No borehole velocity constraints applied. Max velocity = 2819 rn/s; mm velocity = 1724 rn/s. (b) Velocity constraints imposed at boreholes. Max velocity 2801 m/s; mm velocity = 1724 rn/s. 0 0 0 too a 250 0 0 rhorizontal position (m) Fig. 8.6. Reconstructed velocity model (contour interval = 100 m/s) obtained by inverting noise contaminated VSP and crosswell traveltimes. Max velocity = 2810 m/s; nun velocity = 1731 rn/s. 188 Figure 8.6 illustrates that the tornographic inversion is stable when the synthetic tray eltimes are contaminated with small amounts of random noise. Random numbers drawn from a uniform probability distribution on Β±4 ms are added to the exact traveltimes. The iterative inversion is initiated with the same constant slowness model used for the previ ous example. The convergence criterion for terminating iterations is arbitrarily selected to be 2.5 ms of rms traveltime error, or approximately one standard deviation of the noise. No borehole constraints are imposed. An accurate velocity reconstruction results when the weights of the flattening constraints (..jij and in equation (8.14)) are set sufficiently high. C 0 0 250 500 189 The final example examines the ability of the constrained inversion algorithm to recon struct the interwell velocity field when only crosshole traveltimes are available. Figure 8.7 displays the 81 crosswell raypaths. The zone of reduced ray coverage at shallow depths is much more extensive. Images obtained by inverting error free traveltimes are illustrated in Figure 8.8. Once again, both horizontal and vertical flattening constraints are applied. The image in Figure 8.8a displays a shallow low velocity zone and a deeper, dipping high velocity zone. However, the peak of the high velocity anomaly is shifted downdip by a significant distance. An inversion including borehole velocity constraints is displayed in Figure 8.8b. - Surprisingly, the reconstruction is not improved; rather, a spurious region with high velocity gradients develops near the right well. The tomographic inversions in this section are all initiated with uniform slowness models calculated via the method described in Appendix H. Similar results are obtained from a range of nearby slowness values. Alternately, a nonuniform starting-mode1 obtained by horizontal interpolation of the borehole slowness values can be be used. For the examples considered here, this latter technique does not yield any obvious improvement in the velocity reconstructions. 8.8 Conclusion Finite-difference traveltime computation offers an attractive alternative to conventional raytracing for tomographic inversion purposes. The method is sufficiently rapid and accu rate, and handles all of the various wave types that constitute first arrivals. Moreover, since traveltimes are computed throughout a slowness model, very general recording geometries are easily accommodated. The main limitation of the technique is that it is restricted to first arriving waves. Hence, the present formulation cannot be applied to the reflection tomogra phy problem. However, current efforts to generalize finite-difference computation methods to reflection traveltimes (e.g., Podvin and Lecomte, 1991) are encouraging. The introduction of constraint information into traveltime tomography is a responsible way to address the nonuniqueness inherent in this inverse problem. Constraining information may arise from known geological or geophysical properties of the subsurface velocity model D o q I- β. o o de pt h (rn ) 0 N 0 0 0 tI) I-β . 0 C .. o C, β 0 00 . 0 0 Ct, 50 0 25 0 0 0 191 500 Fig. 8.8. Reconstructed velocity models (contour interval = 100 m/s) obtained by inverting only the crosswell traveltimes. (a) No borehole velocity constraints applied. Max velocity = 2811 m/s; mm velocity = 1792 rn/s. (b) Velocity constraints imposed at boreholes. Max velocity = 2862 m/s; mm velocity = 1728 rn/s. 0 - 0 0 0 500250 0 0 tO0 250 horizontal position (m) 192 (say, from outcrops or borehole logs). Alternately, constraints may derive from a desire to impose certain reasonable attributes, like flatness or smoothness, on the constructed model. The method of incorporating constraints into the mathematical inversion procedure is adaptable to either viewpoint. Linear equality constraints are applied directly to the constructed model, rather than to a model perturbation, and are satisfied in the least squares sense. The examples illustrate the imposition of flattening constraints and known borehole information in the reconstruction of a smoothly varying interwell velocity field. Inclusion of this constraint information allows the solution of a problem that is strongly underdetermined - (243 data and 10000 unknowns). Nevertheless, a superior result is not necessarily achieved by the addition of the borehole velocity constraints. 193 CHAPTER 9 SUMMARY This thesis contributes to knowledge in three specific areas of seismic refraction travel- time analysis: (1) Forward modeling and inversion of head wave arrival times is extended to three dimen sional layered earth models in Chapters 2, 3, and 4. This class of models is characterized by uniform velocity layers bounded by plane interfaces with arbitrary strike and dip. Com putational procedures are developed for both single-layer and multilayered models. In the single-layer case, closed form mathematical solutions to the forward and inverse problems exist. A rigorous derivation of a traveltime equation for critically refracted waves propa gating within a multilayered model is also given. The resulting expression is not strictly βclosed formβ; it requires a minimal amount of raytracing to evaluate numerically. This formula forms the basis of an iterative head wave traveltime inversion algorithm designed to recover the parameters defining a single-layer or multilayered earth model. Inclusion of constraint information in the procedure, in the form of inequality relations satisfied by the model parameters, often governs the ability of the algorithm to converge to a realistic solu tion. Tests with simulated and field data, acquired in various recording geometries, indicate that the single-layer version of the algorithm is reasonably robust. However, the multilay ered inversion algorithm appears to require fairly restrictive constraints in order to operate effectively. (2) Improvements to various two dimensional refraction traveltime inversion methods that allow for nonplane interfaces and/or variable velocity media are developed in Chapters 5, 6, and 7. Two existing techniques (the generalized reciprocal method and the wavefront method) are extensively analyzed and a new method (critical offset refraction profiling) is proposed. This new technique is an improvement over the GRM in that it explicitly incorporates undulating interfaces and horizontally varying refractor velocity into the model. Moreover, point values of interface depth, interface dip, and refractor velocity are obtained from the observed arrival times. Hence, a depth profile of the critically refracting horizon can 194 be constructed via interpolation techniques. A rational procedure for calculating a smooth depth profile, based on methods of linear inverse theory, is presented. This latter technique is not limited for use solely with the critical offset refraction profiling method, but may augment any traveltime analysis procedure that yields point depth or dip estimates of the subsurface interface (such as the modified GRM developed here). Finally, an automated implementation of the classical wavefront method for interpreting refraction arrival times indicates that an image of an interface can be obtained directly from the picked times, with no intermediate computational steps. A finite-difference propagation algorithm is used - to downward continue the observed times through a heterogeneous near-surface velocity structure. A simple imaging condition involving the reciprocal time then defines the locus of the subsurface refracting horizon. Tests with synthetic data indicate that both antidinal and synclinal structures can be imaged accurately. Shallow refraction data acquired at an archeological site is also used to assess the workability of the algorithm. (3) Application of finite-difference traveltime computation methods to the two dimensional tomographic inverse problem is treated in Chapter 8. An iterative algorithm is developed for reconstructing a. P-wave velocity field from measured first arrival times. Rapid and accurate forward modeling of all first arrival types (direct waves, head waves, and diffractions) for arbitrary source-receiver geometries is achieved with a finite-difference algorithm. Curved raypaths, needed for converting traveltime residuals into localized updates to the slowness model, are generated by following the steepest descent direction through the computed tray eltime field from each receiver back to the source. Incorporation of constraint information into the procedure, in the form of horizontal and vertical first difference regularization, serves to stabilize the inversion and drive the solution toward a model with a. preferred character (i.e., a βflatβ slowness model). Tests with simulated vertical seismic profile and crosswell arrival times indicate that the algorithm can successfully reconstruct a smoothly varying interwell velocity field, even when the problem is severely underdetermined. These contributions can be applied to the solution of practical problems in seismic refraction exploration, for various objectives and at many different scales. The stage is 195 now set for generalizing the techniques to more complicated three dimensional models with nonpiane interfaces and/or nonuniform velocities. 196 REFERENCES Ackermanu, H.D., Pankratz, L.W., and Dansereau, D., 1986, Resolution of ambiguities of seismic refraction traveltime curves: Geophysics, 51, 223-235. Adachi, R., 1954, On a proof of fundamental formula concerning refraction method of geo physical prospecting and some remarks: Kumamoto J. Sd., Ser. A, 2, 18-23. Aidridge, D.F., Dosso, S.E., Endres, A.L., and Oldenburg, D.W., 1991, New methods for constructing flattest and smoothest models: Inverse Problems, βT, 499-513. Ali Ak, M., 1990, An analytical raypath approach. to the refraction wavefront method: Geophys. Prosp., 38, 971-982. -Barton, D.C., 1929, The seismic method of mapping geologic structure, in Geophysical prospecting: Am. Inst. Mm. Metallurg. Eng., 572-624. Berni, A.J., and Roever, W.L., 1989, Field array performance: theoretical study of spatially correlated variations in amplitude coupling and static shift and case study in the Paris Basin: Geophysics, 54, 451-459. Berry, M.J,, and West, G.F., 1966, An interpretation of the first-arrival data of the Lake Superior experiment by the time-term method: Bull., Seis. Soc. Am., 56, 141-171. Bishop, T.N., Bube, K.P., Cutler, R.T., Langan, R.T., Love, P.L., Resnick, J.R., Shuey, R.T., Spindler, D.A., and Wyld, 11.W., 1985, Tomographic determination of velocity and depth in laterally varying media: Geophysics, 50, 903-923. Bois, P., La Porte, M., Lavergne, M., and Thomas, G., 1972, Well-to-well seismic measure ments: Geophysics, 37, 471-480. Braile, L.W., 1973, Inversion of crustal seismic refraction and reflection data: J. Geophys. Res., 78, 7738-7744. Bregman, N.D., Bailey, R.C., and Chapman, 0.11., 1989, Crosshole seismic tomography: Geophysics, 54, 200-215. Briickl, E., 1987, The interpretation of traveltime fields in refraction seismology: Geophys. Prosp., 35, 973-992. Oerven, V., and Ravindra, R., 1971, Theory of seismic head waves: Univ. of Toronto Press. Chander, R., 1977a, On tracing seismic rays with specified end points in layers of constant velocity and plane interfaces: Geophys. Prosp., 25, 120-124. Chander, R., 1977b, Numerical tracing of critically refracted rays: Geophysics, 42, 1048- 1052. Chin, S.K.L., Kanasewich, E.R., and Phadke, S., 1986, Three-dimensional determination of structure and velocity by seismic tomography: Geophysics, 51, 1559-1571. 197 Chiu, S.K.L., and Stewart, R.R., 1987, Tomographic determination of three-dimensional seismic velocity structure using well logs, vertical seismic profiles, and surface seismic data: Geophysics, 52, 1085-1098. Clayton, R.W., and McMechan, G.A., 1981, Inversion of refraction data by wave field con tinuation: Geophysics, 46, 860-868. V Cunningham, A.B., 1974, Refraction data from single-ended refraction proffles: Geophysics, 39, 292-301. De Amorim, W.N., Hubral, P., and Tygel, M., 1987, Computing field statics with the help of seismic tomography: Geophys. Prosp., 35, 907-919. Diebold, J.B., 1987, Three-dimensional traveltime equation for dipping layers: Geophysics, 52, 1492-1500. Diebold, J.B., and Stoffa, P.L., 1981, The traveltime equation, tau-p mapping, and inversion of common midpoint data: Geophysics, 46, 238-254. Dines, K.A., and Lytle, R.J., 1979, Computerized geophysical tomography: Proc. IEEE, 67, 1065-1073. Dix, C.H., 1935, Note on the theory of seismic prospecting: J. Soc. Petr. Ceophys., 6, 34-43, reprinted in 1947, Early geophysical papers of the Society of Exploration Geophysicists: Soc. Expi. Geophys., 278-287. Dix, C.H., 1941, Notes on refraction prospecting: Geophysics, 6, 378-396. Dooley, J.C., 1952, Calculation of depth and dip of several layers by refraction seismic method, in Thyer, R.F., and Vale, K.R., Geophysical surveys, Oaklands-Coorabin coal- field, New South Wales: Bull., Bur. Mi Res. Austral., Appendix 1. Dunkin, J.W., and Levin, F.K., 1971, Isochrons for a three-dimensional seismic system: Geophysics, 36, 1099-1137. Ewing, M., Wooflard, G.P., and Vine, A.C., 1939, Geophysical investigations in the emerged and submerged Atlantic coastal plain, part III: Barnegat Bay, New Jersey, section: Bull., Geol. Soc. Am., 50, 257-296. Gardner, L.W., 1939, An areal plan of mapping subsurface structure by refraction shooting: Geophysics, 4, 247-259. Gersztenkorn, A., and Scales, J.A., 1988, Smoothing seismic tomograms with alpha-trimmed means: Geophys. J., 92, 67-72. GjΓΈystdal, H., and Ursin, B., 1981, Inversion of reflection times in three dimensions: Geo physics, 46, 972-983. Hadjidaki, E., 1988, Preliminary report of excavations at the harbor of Phalasarna in West Crete: Am. J. Archaeol., 92, 463-479. Hagedoorn, J.G., 1959, The plus-minus method of interpreting seismic refraction sections: Geophys. Prosp., 7, 158-182. 198 Hales, F.W., 1958, An accurate graphical method for interpreting seismic refraction lines: Geophys. Prosp., 6, 285-294. Hampson, D., and Russell, B., 1984, First-break interpretation using generalized linear in version: J. Can. Soc. Expi. Geophys., 20, 40-54. HasselstrΓΆm, B., 1969, Water prospecting and rock-investigation by the seismic refraction method: Geoexploration, 7, 113-132. Hatherly, P.J., 1980, Computer processing of seismic refraction data: Bull., Austral. Soc. Expi. Geophys., 11, 69-74. Hatherly, P.J., 1982, A computer method for determining seismic first arrival times: Geo physics, 47, 1431-1436. llatlierly, P.J., and Neville, M.J., 1986, Experience with the generalized reciprocal method of seismic refraction interpretation for shallow engineering site investigations: Geophysics, 51, 255-265. Hawkins, L.V., 1961, The reciprocal method of routine shallow seismic refraction investiga tions: Geophysics, 26, 806-819. Heiland, C.A., 1929, Modern instruments and methods of seismic prospecting, in Geophys ical prospecting: Am. Inst. Mm. Metallurg. Eng., 625-643. Heiland, C.A., 1940, Geophysical exploration: Prentice-Hall Inc. Hill, N.R., 1987, Downward continuation of refracted arrivals to determine shallow structure: Geophysics, 52, 1188-1198. Hubral, P., 1976, Interval velocities from surface measurements in the three-dimensional plane layer case: Geophysics, 41, 233-242. Hunter, J.A., and Pullan, S.E., 1990, A vertical array method for shallow seismic refraction surveying of the sea floor: Geophysics, 55, 92-96. Johnson, S.H., 1976, Interpretation of split-spread refraction data in terms of plane dipping layers: Geophysics, 41, 418-424. Johnson, L.E., and Gilbert, F., 1972, Inversion and inference for teleseismic ray data, in Bolt, B.A., Ed., Methods of computational physics, 12: Academic Press, 231-266. Jones, G.M., and Jovanovich, D.B., 1985, A ray inversion method for refraction analysis: Geophysics, 50, 1701-1720. Kanasewich, E.R., and Chin, S.K.L., 1985, Least-squares inversion of spatial seismic refrac tion data: Bull., Seis. Soc. Am., 75, 865-880. Kilty, K.T., Norris, R.A., McLamore, W.R., Hennon, K.P., and Euge, K., 1986, Seismic refraction at Horse Mesa dam: an application of the generalized reciprocal method: Geophysics, 51, 266-275. Lankston, R.W., 1989, The seismic refraction method: a viable tool for mapping shallow targets into the 1990s: Geophysics, 54, 1513-1542. 199 Lankston, R.W., 1990, High-resolution refraction seismic data. acquisition and interpretation, in Ward, S.H., Ed., Geotecimical and environmental geophysics, volume 1: review and tutorial: Soc. Expi. Geophys., 45-73. Lankston, R.W., and Lankston, M.M., 1986, Obtaining multilayer reciprocal times through phantoniing: Geophysics, 51, 45-49. Lin, T., 1989, Prestack traveltime inversion for three-dimensional structure: Geophysics, 54, 359-367. Lytle, L.R., and Dines, K.A., 1980, Iterative ray tracing between boreholes for underground image reconstruction: IEEE Trans. Geosci. Remote Sensing, 18, 234-240. Macrides, C.G., Kanasewich, E.R., and Bliaratha, S., 1988, Multiborehole seismic imaging in steam injection heavy oil recovery projects: Geophysics, 53, 65-75. McMechan, G.A., Harris, J.M., and Anderson, L.M., 1987, Cross-hole tomography for strongly variable media with applications to scale model data: Bull., Seis. Soc. Am., 77, 1945-1960. Merrick, N.P., Odins, J.A., and Greenhaigh, S.A., 1978, A blind zone solution to the problem of hidden layers within a sequence of horizontal or dipping refractors: Geophys. Prosp., 26, 703-721. Mota, L., 1954, Determination of dips and depths of geological layers by the seismic refraction method: Geophysics, 19, 242-254. Nolet, G., 1985, Solving or resolving inadequate and noisy tomographic systems: J. Comp. Phys., 61, 463-482. Nolet, G., 1987, Seismic wave propagation and seismic tomography, in Nolet, G., Ed., Seismic tomography with applications in global seismology and exploration geophysics: D. Reidel Pubi. Co., 1-23. Ocola, L.C., 1972, A nonlinear least-squares method for seismic refraction mapping - part II: model studies and performance of refra.map method: Geophysics, 37, 273-287. Oldenburg, D.W., 1981, A comprehensive solution to the linear deconvolution problem: Geophys. J. Roy. Astr. Soc., 65, 331-357. Oldenburg, D.W., 1984, An introduction to linear inverse theory: IEEE Trans. Geosci. Remote Sensing, 22, 665-674. Olsen, K.B., 1989, A stable and flexible procedure for the inverse modelling of seismic first arrivals: Geophys. Prosp., 37, 455-465. Paige, C.C., and Saunders, M.A., 1982, LSQR: an algorithm for sparse linear equations and sparse least squares: ACM Trans. Math. Software, 8, 43-71. Palmer, D., 1980, The generalized reciprocal method of seismic refraction interpretation: Soc. Expi. Geophys. 200 Palmer, D., 1981, An introduction to the generalized reciprocal method of seismic refraction interpretation: Geophysics, 46, 1508-1518. Palmer, D., 1986, Refraction seismics, the lateral resolution of structure and seismic velocity: Geophysical Press. Palmer, D., 1990, The generalized reciprocal method - an integrated approach to shallow refraction seismology: Expi. Geophys., 21, 33-44. Palmer, D., 1991, The resolution of narrow low-velocity zones with the generalized reciprocal method: Geophys. Prosp., 39, 1031-1060. Parker, R.L, Shure, L., and Hildebrand, J.A., 1987, The application of inverse theory to seamount magnetism: Rev. Geophys., 25, 17-40. Pavlenkin, A.D., Saiculina, S., and Vinnik, A.A., 1986, Method of conjugate points in the interpretation of traveltime curves of refracted waves: Izvestia, Earth Physics, 22, 795- 803. Phadke, S., and Kanasewich, E.R., 1990, Seismic tomography to obtain velocity gradients and three-dimensional structure and its application to reflection data on Vancouver Island: Can. J. Earth Sci., 27, 104-116. Podvin, P., and Lecomte, I., 1991, Finite difference computation of traveltimes in very con trasted velocity models: a massively parallel approach and its associated tools: Geophys. J. Tnt., 105, 271-284. Press, W.H., Flannery, B.P., Teukolsky,S.A., and Vetterling, W.T., 1986, Numerical recipes, the art of scientific computing: Cambridge Univ. Press. Radcliffe, R.D., Balanis, C.A., and Hill, Jr., H.W., 1984, A stable geotomography technique for refractive media: IEEE Trans. Geosci. Remote Sensing, 22, 689-703. Reiter, L., 1970, An investigation into the time term method in refraction seismology: Bull., Seis. Soc. Am., 60, 1-13. Richards, P.G., Witte, D.C., and Ekstr6m, G., 1991, Generalized ray theory for seismic waves in structures with planar nonparallel interfaces: Bull., Seis. Soc. Am., 81, 1309-1331. Rockwell, D.W., 1967, A general wavefront method, in Musgrave, A.W., Ed., Seismic re fraction prospecting: Soc. Expi. Geophys., 363-415. Russell, D.R., Keller, G.R., and Braile, L.W., 1982, A technique to determine the three- dimensional attitude and true velocity of a refractor: Geophysics, 47, 1331-1334. Scales, J.A., 1987, Tomographic inversion via the conjugate gradient method: Geophysics, 52, 179-185. Scales, J.A., Docherty, P., and Gersztenkorn, A., 1990, Regularisation of nonlinear inverse problems: imaging the near-surface weathering layer: Inverse Problems, 6, 115-131. Scheidegger, A.E., and Willmore, P.L., 1957, The use of a least squares method for the interpretation of data from seismic surveys: Geophysics, 22, 9-21. 201 Schenck, F.L., 1967, Refraction solutions and wavefront targeting, in Musgrave, A.W., Ed., Seismic refraction prospecting: Soc. Expi. Geophys., 416-425. Schwarz, S.D., 1990, Detection of destressed rock and potential collapse features above old mine workings by the seismic refraction method, in Ward., S.H., Ed., Geotechnical and environmental geophysics, volume 3: geotechnical: Soc. Expi. Geophys., 281-287. Shah, P.M., 1973, Ray tracing in three dimensions: Geophysics, 38, 600-604. SjΓΆgren, B., 1980, The law of parallelism in refraction shooting: Geophys. Prosp., 28, 716-743. Smith, T.J., Steinhart, J.S., and Aldrich, L.T., 1966, Lake Superior crustal structure: J. Geophys. Res., 71, 1141-1172. -Sorrells, G.G., Crowley, J.B., and Veith, K.F., 1971, Methods for computing ray paths in complex geological structures: Bull., Seis. Soc. Am., 61, 27-53. Tarrant, L.H., 1956, A rapid method of determining the form of a seismic refractor from line profile results: Geophys. Prosp., 4, 131-139. Thornburgh, H.R., 1930, Wave-front diagrams in seismic interpretation: AAPG Bull., 14, 185-200. Vidale, J., 1988, Finite-difference calculation of travel times: Bull., Seis. Soc. Am., 78, 2062-2076. Vidale, J.E., 1990, Finite-difference calculation of traveltimes in three dimensions: Geo physics, 55, 521-526. White, D.J., 1989, Two-dimensional seismic refraction tomography: Geophys. J., 97, 223- 245. Wiffiamson, J.H., 1968, Least-squares fitting of a straight line: Can. J. Phys., 46, 1845-1847. Willmore, P.L., and Bancroft, A.M., 1960, The time term approach to refraction seismology: Geophys. J. Roy. Astr. Soc., 3, 419-432. Wyrobek, S.M., 1956, Application of delay and intercept times in the interpretation of multilayer refraction time distance curves: Geophys. Prosp., 4, 112-130. York, D., 1966, Least-squares fitting of a straight line: Can. J. Phys., 44, 1079-1086. York, D., 1967, The best isochron: Earth Plan. Sci. Lett., 2, 479-482. York, D., 1969, Least squares fitting of a straight line with correlated errors: Earth Plan. Sci. Lett., 5, 320-324. Zelt, C.A., 1989, Seismic structure of the crust and upper mantle in the Peace River Arch region: Ph.D. thesis, Univ. of British Columbia. Zelt, B.C., Effis, R.M., and Clowes, R.M., 1990, SCoRE β89; the Southern Cordillera re fraction experiment; description of data set: Lithoprobe report, 20, Univ. of British Columbia. 202 APPENDICES Appendix A: Conic sections in poiar coordinates Expressions for conic sections in plane polar coordinates can be found in most textbooks on analytic geometry. However, these formulae typically assume that the coordinate origin is located either at the center of symmetry or at a focus of the curve. In this appendix, an alternate representation of the effipse or hyperbola in polar coordinates is derived. The principle axes of the conic are rotated by an arbitrary angle S with respect to the coordinate axes, and the origin resides at any point along the major axis of symmetry. In Figure Al, the primed coordinate system is obtained by a clockwise rotation of the unprimed system through the angle 9. The coordinates of a given point relative to each system are related via = x cos 9 + y sinS, yβ = βz sin9 + y cos 9. (Al) Consider an ellipse with center located on the zβ axis at (xβ, yβ) = (la, 0). The sexuimajor axis a and the semiminor axis b are parallel to the xβ and yβ coordinate axes, respectively. Figure Al depicts the case where a < l. However, the situation a> 1C is also allowed; in the case, the effipse encompasses the origin. In the primed reference frame, the equation of the effiptical locus is (z-lc)2 + (c)2 = i. (A2) Substituting from relations (Al) yields (x cos s + y sinS β ic)2 + ( cos x sin 9)2 = 1. 203 However b2 = a2(l β e2) where e is the ellipse eccentricity (0 < e < 1). Thus (1βe)(xcos6+ysin6βl) + (ycos8βxsin9)2 = a2(1βe). This expression is converted to plane polar coordinates (X, a) via the substitutions x X cos a and y = X sin a. After some algebraic reduction, a quadratic form in the radial coordinate X is obtained: X2 [i β e2 cos2(a β 0)] β 2X [(1 β e2)l cos(a β 6)] + [(1 βe2)(l β a2)j = 0. (A3) The quadratic formula can be used to solve equation (A3) for X as a function of the azimuthal angle a. However, this is algebraically tedious as well as unecessary for present purposes. Note that the conventional expressions for an ellipse in polar coordinates are obtained as special cases of equation (A3). Thus, setting l = 0 (i.e., origin coincident with the ellipse center) immediately yields I 1βe2 X(a) = a41 V lβecos(aβ0) Also, if l = a e, then the origin is located at a focus of the ellipse. In this case (A3) gives 1βe2 X(a) = This form is common in mechanics; it describes the trajectory of a particle moving under the influence of an inverse square central force. A similiar analysis reveals that equation (A3) also applies to a hyperbola with the same center location and principle axis orientation. The only difference is that the numerical value of the eccentricity exceeds one (e > 1) for a hyperbola. 204 x x, 0 y y, Fig. Al. A noncentered, rotated effipse. 205 Appendix B: Traveltime analysis for a simple 3D model A compact formula for head wave traveltime in a simple three dimensional earth model is derived in Chapter 2 (equation (2.10)). The purpose of this appendix is to verify that the more general expression (3.12) reduces to this known result in the special case of a two media model with a horizontal surface. Evaluating equation (3.12) with 1,z = 0, mi,z 1, and with the interface index k = 2 gives - T2(xs, , ) [cos 1I q12,z + sin β q12,y] x + hi(O, O)(p12,z β q12,z) β xS(p12,zβ q12,z) β ys(p12,i,β q12,) (B1) Expressions for the zyz components of the unit propagation vectors are now determined. As indicated in Chapter 2, the entire head wave raypath for this situation (depicted in Figure 2.2) is confined to a single plane. The two vectors n2 and p22 form an orthonormal basis for all vectors in this plane. Hence, the propagation vectors p12 and q12 can be resolved along this basis. Geometric analysis of Figure 2.2 yields P12 (cos i)n2 + (sin ic)P22, qrz (β cos i)n2 + (sin ic)p22, (B2) where i is the critical angle. The unit normal n2 to the refracting interface is given by = (sin 4 cos 02)i + (sin 4 sin02)J + (cos qfj)k. (B3) An expression for the critically refracted propagation direction vector P22 is obtained by recognizing that rQ β rp P22 = trQβrp1I 206 where rp and r are the position vectors of points P and Q in Figure 2.2. Thus rq β rp = (rR + dRn2) β (rs + dsn2) = (XRβZS)I+(YRβYS)i+(dRβdS)n2 = (X cos Pβ)i + (X sin [β)j + (βX sin 8)112, where equation (2.7) has been used and the angle S is defined by sin S = sin 4i cos(βJβ β 82). Also, equation (2.9) implies rQ β rp = X cos S. Hence (cos P)i + (sin βIβ)j β (sin S)n2 P22 = . (B4)cos Substituting these expressions for 2 and P22 into equations (B2) yields the rather formidable formulae sin #2 COS 82 cos(i + 8) + sin i cos P12 1cos S + sin #2 sin 82 cos(i + 8) + sinicsinβ11 + [cos #2 cos(i + 8) k, (B5)cosS j L cosS and β Sfl #2 COS 82 cos(i 8) + Sin ic COS iβ 1. q12= coso J + βsinq52sinO2cos(iCβ8)+sinisin1 . + cos#2os(i8) k (B6) cos8 cosS It is straightforward to verify that II P12 11=11 qi 11= 1, P12 2 = cosi, and q β’n2 = β cos i, as is required. Moreover, equations (B5) and (B6) imply that 207 P12,z β q12,z 2 COS i Sifl #2 COS 62, P12,,, q12,y = 2 COS ic Sifl #2 sin 02, P12,z β q12,z 2 cos cos #2, cos E + sin β1 q12,y = sin(i β 8). These relations are substituted into the traveltime expression (B1) to obtain the final result T2(ZS,YS,X, IIβ) = sin(icβ + 2dscosi (B7) where ds is the perpendicular distance from the source S to the refracting interface: ds =h1(0,O)cos β sin(zScosO2 +yssin62). Except for some minor notational changes, equation (B7) agrees exactly with the traveltime formula (2.10) previously developed for this simple model situation. 208 Appendix C: Traveltime analysis for a simple 2D model A useful check on the validity of the general head wave traveltime formulae in Chapter 3 is provided by analyzing a specific situation for which a closed form traveltime solution exists. This is especially important because several of these expressions disagree with analogous formulae published recently by Diebold (1987). The inconsistency between the two results is evident when source and receiver are located on separate interfaces with different dip angles. In this case, it is not possible to effect a coordinate frame rotation such that both interfaces become horizontal (i.e., parallel to the zy plane). Consider the simple two dimensional earth model depicted in Figure Cl. Two subsurface interfaces have the same dip angle Γ§. The source S is located at the coordinate origin on the surface and the receiver R is located on subsurface interface 2. L is the source-receiver range measured parallel to the dipping interfaces. The perpendicular distances from the origin to the subsurface interfaces are d2 and d3; vertical layer thicknesses at the same point are h1 and h2. The traveltime of a head wave formed on interface 3 can be derived from first principles of refraction traveltime analysis. Rotating the perspective through the small angle cp trans forms the situation into a one dimensional problem. The total traveltime is easily obtained by summing the refraction delay time associated with each layer that the wave traverses, together with L/v3. The result is T3 = . + d2cos9j3 + 2(d3β)cos93 (Cl) where sin 9jj = v/v1. The following geometric relationships are evident from Figure Cl: d2 = h1 cos y, d3 β d2 = h2 cos , L = + d2 tan Γ§a, cos Γ§o 209 / Fig. Cl. A simple 2D earth model with three media. S is a surface source and R is a buried receiver. where ZR is the horizontal coordinate of the receiver. Substituting these expressions into equation (Cl) yields T3(ZR) = ZR + h1 cos(013 β + 2h cos 823 COS V3COS9 1)1 V2 (C2) The traveltime predicted by the general formula (3.16) is now compared with this specific traveltime solution. For the situation being examined, (3.16) simplifies considerably. The following conditions hold: S d1 XR / d2 vi I / / R V3 (p 210 (i) XS=YS=YR=O. (ii) j = 1, 1 = 2, k = 3. (iii) h1(O,0) = h1, h2(O,O) = h2. (iv) = β Sifl 2,z = COS 50. Evaluating formula (3.16) with these parameters yields T3(R) = XR COS 50 q23, + SIR 50 q23,z + h1 P13,z + h2(p3, β q23,z) (C3) V2COSSO V2 Since the earth model is strictly two dimensional, and the recording proffle is oriented normal to the common strike direction of the interfaces, the unit propagation direction vectors are contained entirely within the xz plane. Thus, the components of these vectors can be obtained by further geometric analysis of Figure Cl: P13,z cos(6i β so), P23,z cos(623 β q23,z sin(623 + so), q23,z = β cos(623 + 50). Substituting these expressions into equation (C3) gives ZR Sm 623 h1 cos(613 β ) 2h cos COS 50 T3(ZR)= + + v2cosΓ§o v1 V2 But since sin 923 = v2/v3, this equation reduces immediately to = ZR + h1 cos(913 β 0) + 2h COS 023 COS (C4) vacOsΓ§o vi V2 211 This is in complete agreement with the known solution (C2)! In contrast, equation (21) in Diebold (1987), when applied to the model depicted in Figure Cl, becomes T3(ZR) = ZR sin(8:3 + + h1 cos(913 β + 2h cos::3 CoScp (C5) Clearly, this expression differs from the expected solution (C2). It reduces to the correct traveltime only if the dip angle Γ§o equals zero, i.e., all interfaces are horizontal. 212 Appendix D: Derivation of formula (6.13) A standard technique from the calculus of variations is used to demonstrate the vadidity of equation (6.13). If the function mβ(z) is perturbed by an arbitrary, but small, amount Smβ(x), then a variation S is induced in the functional 4(mβ). An expression for this variation is obtained by evaluating equation (6.12) at mβ + Smβ. Hence 5cJ (mβ + Smβ) β = 2j {1Lw(z)2m(x) β [eob3 β eprd(mβ)jTWTWp(x;cl,c2)}SmI(x) dx + O(Smβ)2. (D1) (mβ) is extremized by the particular mβ(x) such that S is second order in Smβ. Thus, the integral in equation (Dl) must vanish identically. Since Smβ(x) is an arbitrary small perturbation, this can be achieved only if mβ(x) = --- [cobs β erd(mβ)] WTW P(z;cl,:2) (D2) The mβ(z) that extremizes (mβ) is a linear combination of the functions pj(x; Cl, c2)/w(x). Moreover, the coefficient vector in this linear superposition is proportional to the difference between observed and predicted data. Parker et al. (1987) obtain an analogous result for a different inverse problem formulation. 213 Appendix E: Vanishing of two derivatives This appendix proves that the extremum of the objective function in equation (6.17) is invariant with respect to the two abscissae c and C2. Differentiating 4β with respect to these parameters yields the two equations = aβ [2r WTW + I] a ad ad + 2 [ctlh(b) β d2 [k(b) β (b β ci)h(b)j β eobs] a β2dh(b)TW[dlh(b) β d2 {k(b) β (b β ci)h(b)] β e0b3 + Fa] (El) and = aT [2r WTW + u ij c9d2 + 2 [cilh(b) β d2 [k(b) β (b β ci)h(b)] β eobs] T a, (E2) where Iβ stands for F(ci, c2). These expressions are simplified by using the previous equation (6.18) obtained by extreniizing the objective function with respect to the coefficient vector a: Pa + ,L(WTW)_la e03 β dih(b) + d2 [kb β (b β ci)h(b)j. (E3) Substituting (E3) into (El) and (E2) immediately yields the simpler forms = ,taT + 2,udh(b)Tcr, (E4) 214 Tβ C2 C2 Expressions for the partial derivatives of the inner product matrix F are required. Dif ferentiating equation (6.16) with respect to Cl yields = L w(x)_2 pT + p (9P)T] dx, (E6) where p stands for p(z; c,c2). An analogous formula is obtained for bF/8c2.The derivatives of the kernel function vector p(x; ci, c2) are obtained by differentiating expression (6.10). Hence c2) = h(b) [H(x - ci) - - c2)j + h(b) (x - Cl) S(x - Cl), (E7) OP(x;::,c2) = [k(b) β (b β ci)h(b)jS(z β c2). (E8) Substituting (E7) into (E6) and integrating gives the expression = h(b)qT + qh(b)T, (E9) where q is an auxilliary vector defined by q f2w(z)p(x;c, c2) dx. Similarly, substi tuting (E8) into the analogue of (E6) for ΓF/0c2 gives the result 0C2 = w(c2)2 [[k(b) β (bβ ci)h(b)]p(c2;ci,c)T + p(c2;cl,c2) {k(b) β (b β cl)h(b)]T]. (E10) 215 Relations (E9) and (Elo) are the required formulae. Substituting these expressions into (E4) and (E5) and reducing yield = 2p (d2 β aTq) [aTh(b)], (Eli) 8C2 = β2p w(c2 )_2 [a Tp(c; c,c2)] [a Tk(b) β (b β ci) a Th(b)] (E12) However, if the objective function 4 is already extremized with respect to the two parameters d1 and d2, then aTh(b) = aTk(b) = 0. Thus, equations (Eli) and (E12) reduce to the desired results 0, β 0. (E13) Note that both of the conditions a Th(b) = 0 and a Tk(b) = 0 are necessary for 8c1/8 2 to vanish. 216 Appendix F: Spatially correlated traveltime errors Suppose that x is a vector of ii random variables with covariance matrix C. Then, the vector of n random variables y given by the linear transformation y = Ax possesses the covariance matrix C1, = ACXAT. (Fl) If C,, and C1, are prescribed, then equation (Fl) can be solved for the (n x n) transformation matrix A. Since the covariance matrices are symmetric and positive definite, they may be factored via Cholesky decomposition as follows: i _y r β .LJ 12Z β-β1, β U1, .L J1,, where L,, and L1, are lower triangular matrices. Then, one solution of (Fl) is A = L1,;, (F2) which can be readily verified by substitution. An important special case occurs when the z1 are a set of independent random variables drawn from the one dimensional normal distribution with zero mean and unit standard deviation. In this situation, the covariance matrix C,, equals the identity matrix. The transformed random variables are given simply by y = L1, x. Moreover, each y is normafly distributed because it is a linear combination of independent normal variates. Within the context of the arrival time picking problem, each y is considered to be a r.ndom time error with an assigned standard deviation o. A double-tailed exponential function is used to evaluate the correlation coefficient between the arrival time errors at 217 geophone stations i and j (Berni and Roever, 1989). The elements of the covariance matrix become {cgj = exp(βd/D)ojoj, (F3) where djj is the distance between stations i and j, and D is an adjustable parameter called the correlation distance. A large value for D (relative to the dj) implies that the individual arrival time errors are highly correlated, and vice versa. Standard FORTRAN subroutines (Press et aL, 1986, p. 192-199) are used to generate independent normal variates aj with zero mean and unit variance. Hence, the y calculated via the above procedure become a set of correlated, normally distributed traveltime errors. 2i8 Appendix G: Ray initiation and termination logic If a receiver is located on a grid node, on a horizontal grid line, or on a vertical grid line (points 1, 2, and 3 in Figure Cia, respectively), then the average of the steepest descent angles from the immediately adjacent cells determines which cell the raypath enters first. However, if the receiver is located within a cell, as in Figure Cib, then the logic contained in Table Ci is used to determine the initial path segment and the cell that the raypath subsequently enters. Angles /9i through /38 are all positive acute angles and are illustrated in Figure Gib. Now assume that points 1, 2, and 3 are sources, rather than receivers. Then, the set of immediately adjacent cells constitutes a near source zone. Tracing of an incoming ray continues until the boundary of this zone is encountered. For example, the four cells surrounding point 1 in Figure Cia form a square with side length 2h. For source points 2 and 3, the near source zones are rectangles with vertical and horizonal dimensions of 2h x h and h x 2h, respectively. Finally, if the source is located within a cell (Figure Gib), then iterative ray tracing proceeds to the boundary of the enclosing square cell with side length h. In all cases, the final raypath segment is obtained by drawing a straight line from the boundary point directly to the source position. Zi h Zr Zj+1 3 a h 1 .2 h h h 219 Fig. Gi (a) Points 1, 2, and 3 denote receivers (or sources) located on a grid node or grid lines. (b) A receiver located within a square grid cell. X1 Xr xi+1 h 220 Angle Range Exit Coordinates Next Cell Indices o <Bk <13i z,. + β x)tan 8k i + 1 j = f3 i + 1 j + 1 - ir/2β/2<Ok<lr/2+1 Zr+(zj+1zr)COtOk Z3f1 i 3+1 Bk=lr/+/3 Zj+1 iβi j+1 1β/4<Bk<+/5 Z zT+(zjβxf)tanOk iβi j Sk=1+/35 Z Zj iβi jβ1 31r/2β/6<8k <31r/2+/37 xr+(zjβz,.)cotOk zi I jβ1 Sk=31r/2+/ i+1 jβi 2ir β /38 <Bk <2w xj+1 z,. + (zjΒ± β zr)tan 0k i + 1 j Table Gi. Ray initiation logic for a receiver located within a cell. Column 1 gives possible ranges for the steepest descent direction 6k assigned to the cell. Angles /9j through /38 are illustrated in Fig. Gib. Columns 2 and 3 give the horizontal and vertical coordinates of the ray exit point on the cell boundary, respectively. (z,., Zr) are the receiver coordinates. Columns 4 and 5 give the indices of the next cell that the ray enters. 221 Appendix H: Optimum starting slowness In the absence of a priori information, a. uniform slowness model is used to initiate the tomographic inversion procedure. The value of this slowness can be chosen to minimize the rms traveltime residual on the first iteration. This tends to reduce the total number of iterations required for convergence to a specified misfit level. Let t be the th observed traveltime, and let the straight line distance between the associated source and receiver be dj. Then, the difference between the observed traveltime and the time predicted by a uniform slowness s is t β s dj. The 12 norm (squared) of all such differences is 4(s) = (td)2 where N is the total number of observed arrival times. Extremizing 4(s) with respect to 8 yields - - EfL1jt sβstβr Alternately, the problem can be posed in terms of an optimum velocity instead of an optimum slowness. The result is v = 1/s. For the combined VSP and crosswell datasets used in Chapter 8, vt equals 2191 rn/s (for both exact and error contaminated traveltimes). For the crossliole data only, v0 = 2310 m/s.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Analysis and inversion of seismic refraction traveltimes
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Analysis and inversion of seismic refraction traveltimes Aldridge, David Franklin 1992
pdf
Page Metadata
Item Metadata
Title | Analysis and inversion of seismic refraction traveltimes |
Creator |
Aldridge, David Franklin |
Date | 1992 |
Date Created | 2008-12-20 |
Date Issued | 2008-12-20 |
Description | Techniques for forward modeling and inversion of head wave traveltimes within the framework of one and two dimensional earth models are well developed. The first portion of this thesis extends these methods to encompass three dimensional layered models. Each critically refracting horizon of the model is approximated by a plane interface with arbitrary strike and dip. An advantage of this simple representation is that rapid computation of head wave traveltimes for arbitrary source-receiver geometries can be achieved with a minimum of ray tracing. Inversion methods are then developed for estimating the parameters defining single-layer and multilayer earth models. For the single-layer model, an algebraic solution to the inverse problem exists if refraction traveltimes are observed along two independent line profiles. For multilayer models and/or nonprofile recording geometries, the inversion is formulated as a constrained parameter optimization problem and solved via linear programming. Inclusion of constraints, in the form of inequality relations satisfied by the model parameters, often governs the ability of the algorithm to converge to a realistic solution. The procedure is tested with traveltimes recorded on broadside profiles in a deep crustal seismic experiment. The second part of this thesis provides specific improvements to various two dimensional refraction traveltime inversion techniques. The generalized reciprocal method (GRM) is reformulated on the basis of an earth model characterized by vertical, rather than normal, layer thicknesses. This allows point values of interface depth to be inferred from the observed traveltimes. A novel interpretation method (critical offset refraction profiling) is described that yields point values of interface depth, interface dip, and refractor velocity. A smooth depth profile of the refracting horizon is then constructed using techniques of linear inverse theory. Finally, an automated version of the classical wavefront method for interpreting refraction traveltimes is developed. Recorded arrival times are downward continued through a near surface heterogeneous velocity structure with a finite-difference propagation algorithm. The locus of a refracting horizon is then obtained by applying a simple imaging condition involving the reciprocal time (the source-to-source traveltime). The method is tested, apparently successfully, on a shallow refraction dataset recorded at an archeological site. The final portion of this thesis develops an iterative tomographic inversion procedure for reconstructing a two dimensional P-wave velocity field from measured first arrival times. Two key features of this technique are (i) use of a finite-difference algorithm for rapid and ac curate forward modeling of traveltimes, and (ii) incorporation of constraint information into the inversion to restrict the nonuniqueness inherent in large scale, nonlinear tomographic inverse problems. Analysis of a simulated vertical seismic profile (VSP) plus crosswell experiment indicates that the inversion algorithm can accurately reconstruct a smoothly varying interwell velocity field. Inclusion of constraints, in the form of horizontal and vertical first-difference regularization, allows the solution of a traveltime tomography problem that is otherwise severely underdetermined. |
Extent | 4060049 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | Eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2008-12-20 |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052957 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of |
Degree Grantor | University of British Columbia |
Graduation Date | 1992-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
URI | http://hdl.handle.net/2429/3268 |
Aggregated Source Repository | DSpace |
Download
- Media
- [if-you-see-this-DO-NOT-CLICK]
- [if-you-see-this-DO-NOT-CLICK]
- ubc_1992_spring_aldridge_david_franklin.pdf [ 3.87MB ]
- [if-you-see-this-DO-NOT-CLICK]
- Metadata
- JSON: 1.0052957.json
- JSON-LD: 1.0052957+ld.json
- RDF/XML (Pretty): 1.0052957.xml
- RDF/JSON: 1.0052957+rdf.json
- Turtle: 1.0052957+rdf-turtle.txt
- N-Triples: 1.0052957+rdf-ntriples.txt
- Original Record: 1.0052957 +original-record.json
- Full Text
- 1.0052957.txt
- Citation
- 1.0052957.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
United States | 20 | 4 |
Norway | 15 | 0 |
France | 8 | 0 |
India | 7 | 0 |
Italy | 7 | 0 |
Sweden | 4 | 0 |
Indonesia | 3 | 0 |
Japan | 3 | 0 |
Egypt | 3 | 2 |
Canada | 3 | 2 |
United Kingdom | 3 | 0 |
Greece | 3 | 0 |
United Arab Emirates | 2 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 54 | 13 |
St Louis | 8 | 0 |
Stockholm | 4 | 0 |
Tokyo | 3 | 0 |
Boulder | 3 | 1 |
Camerino | 3 | 0 |
Dubai | 2 | 0 |
Delhi | 2 | 0 |
Mountain View | 2 | 1 |
Beijing | 2 | 0 |
Cairo | 2 | 1 |
Ashburn | 2 | 0 |
Gaborone | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052957/manifest