ANALYSIS AND INVERSION OF SEISMIC REFRACTION TRAVELTIMESbyDAVID FRANKLIN ALDRIDGEB.A., Amherst College, 1975M.A., The University of California, Berkeley, 1980A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIESDepartment of Geophysics and AstronomyWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAApril 1992© David Franklin Aldridge, 1992Signature(s) removed to protect privacyIn presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)__________________________Department of wc/ A%cQ4o,.17 /The University of British ColumbiaVancouver, CanadaDate 27DE-6 (2/88)Signature(s) removed to protect privacy11ABSTRACTTechniques for forward modeling and inversion of head wave traveltimes within theframework of one and two dimensional earth models are well developed. The first portionof this thesis extends these methods to encompass three dimensional layered models. Eachcritically refracting horizon of the model is approximated by a plane interface with arbitrarystrike and dip. An advantage of this simple representation is that rapid computation of headwave traveltimes for arbitrary source-receiver geometries can be achieved with a minimumof ray tracing. Inversion methods are then developed for estimating the parameters definingsingle-layer and multilayer earth models. For the single-layer model, an algebraic solutionto the inverse problem exists if refraction traveltimes are observed along two independentline profiles. For multilayer models and/or nonproffle recording geometries, the inversion isformulated as a constrained parameter optimization problem and solved via linear programming. Inclusion of constraints, in the form of inequality relations satisfied by the modelparameters, 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 crustalseismic 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 theobserved traveltimes. A novel interpretation method (critical offset refraction profiling) isdescribed that yields point values of interface depth, interface dip, and refractor velocity. Asmooth depth profile of the refracting horizon is then constructed using techniques of linearinverse theory. Finally, an automated version of the classical wavefront method for interpreting refraction traveltimes is developed. Recorded arrival times are downward continuedthrough a near surface heterogeneous velocity structure with a finite-difference propagationalgorithm. The locus of a refracting horizon is then obtained by applying a simple imagingcondition involving the reciprocal time (the source-to-source traveltime). The method isUItested, apparently successfully, on a shallow refraction dataset recorded at an archeologicalsite.The final portion of this thesis develops an iterative tomographic inversion procedurefor 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 accurate forward modeling of traveltimes, and (ii) incorporation of constraint information intothe inversion to restrict the nonuniqueness inherent in large scale, nonlinear tomographicinverse problems. Analysis of a simulated vertical seismic profile (VSP) plus crosswell experiment indicates that the inversion algorithm can accurately reconstruct a smoothly varyinginterwell velocity field. Inclusion of constraints, in the form of horizontal and vertical first-difference regularization, allows the solution of a traveltime tomography problem that isotherwise severely underdetermined.ivTABLE OF CONTENTSAbstract . iiTable of contents ivList of tables . . . . viiiList of figures1 Introduction2 Head wave traveltimes in a three dimensional single-layer earth2.1 Introduction2.2 Earth model and recording geometry2.3 Traveltime derivation2.4 Traveltime inversion2.4.1 Special cases2.4.2 Arbitrary line azimuths2.5 Traveltime isochrons2.6 Inversion of traveltime isochrons2.7 Conclusion3 Head wave traveltimes in a three dimensional multilayered earth3.1 Introduction3.2 Earth model3.3 Raypath geometry3.4 Traveltime derivation3.4.1 Downgoing traveltime3.4.2 Upgoing traveltime3.4.3 Critically refracted traveltime3.4.4 Total traveltime3.4.5 Variants of the basic formula101010• . . 13171820• . . 212627292930• • . 3236• . . 3639• . 40• . 40• . . 423.5 Generalizations of the traveltime formula .3.5.1 Source and receiver on separate interfaces3.5.2 Arbitrary source and receiver locations3.5.3 Arbitrary reference points for layer thickness3.5.4 Reflection traveltime3.6 Rapid traveltime computation3.6.1 General description3.6.2 Raytracing technique3.7 Forward modeling examples3.7.1 Profile geometry3.7.2 VSP geometry3.8 ConclusionV4444464852545557616164664 Inversion of head wave traveltimes for three dimensionalplanar structure4.1 Introduction4.2 Inversion mathematics4.2.1 General theory4.2.2 Calculation of sensitivities .4.3 Synthetic examples4.3.1 Single layer4.3.2 Multiple layers4.4 Field data example4.4.1 Peace River Arch broadside data4.4.2 Static corrections4.4.3 Inversion results4.5 Conclusion686869697276768081818387945 Point depth and dip estimates from refraction traveltimes5.1 Introduction96965.2 Specialization to a 2D model5.3 Derivation of GRM parameters5.3.1 Velocity analysis function5.3.2 Apparent horizontal velocity5.3.3 Generalized time-depth5.3.4 Depth conversion factor5.3.5 Time-depth near a shotpoint5.4 Critical offset refraction profiling5.4.1 Earth model5.4.2 Inversion method5.4.3 Critical offset determination5.5 Conclusionvi971011011031081101121131131151181207 Refractor imaging using an automated wavefrontreconstruction method7.1 Introduction7.2 Finite-difference traveltimes7.2.1 Wavefront construction1221221241251271281301331341371446 Construction of a smooth refractor depth profile6.1 Introduction6.2 Smoothest model construction6.2.1 Modified data equation6.2.2 Objective function6.2.3 Extremizing the objective function . .6.2.4 Constructing the model6.3 Interpolation via the smoothest model6.3.1 General theory6.3.2 Numerical example6.4 Conclusion146146147147vu7.2.2 Wavefront reconstruction7.3 Refractor imaging7.4 Refractor velocity estimation7.5 Field data example7.6 Conclusion1521521601621678 Two dimensional tomographic inversion withfinite-difference traveltimes8.1 Introduction8.2 General theory and model representation8.3 Forward modeling8.4 Ray generation8.5 Inversion mathematics .8.6 Model constraint equations8.7 Synthetic examples8.8 Conclusion1699 SummaryReferencesAppendicesA Conic sections in polar coordinatesB Traveltime analysis for a simple 3D modelC Traveltime analysis for a simple 2D modelD Derivation of formula (6.13)E Vanishing of two derivativesF Spatially correlated traveltime errors .G Ray initiation and termination logic .H Optimum starting slowness193196202202205208212213216218221169171172175178182184189VIIILIST OF TABLES8.1 Example of ray tracing logic 177Cl Ray initiation logic for a receiver located within a cell 220ix2.1 Single layer earth model and coordinate reference frame .2.2 Plan view of surface recording geometry2.3 Head wave raypath2.4 Normalized apparent refractor velocity2.5 Head wave traveltime isochrons3.1 Orientation of i interface of multilayered model- 3.2 Schematic diagram of a raypath critically refracted at kth interface3.3 Snell’s law of refraction at th interface3.4 Snell’s law of reflection at kthl interface3.5 Schematic representation of a critically reflected/refracted raypath3.6 Direct and head wave traveltimes recorded by four reversed profiles3.7 Head wave traveltimes recorded in a VSP configuration4.1 Plan views of two areal recording arrays4.2 Broadside recording geometry for Peace River Arch (PRA) experiment4.3 Receiver static functions for PRA lines A and B4.4 Comparison between PRA observed and predicted traveltimes4.5 Plan view of vertical depths to Moho in PRA region4.6 Vertical depths to Moho beneath PRA lines A and B5.1 Critically refracted raypaths for ORM analysis 1025.2 Comparison of ORM apparent velocities with true velocity 1065.3 Two dimensional model used for critical offset refraction profiling 1145.4 Expanded view of triangle ABC 1175.5 Determination of critical offset distance 1196.1 Refractor elevation model constructed from accurate depth samples 1386.2 Refractor elevation model with optimized parameters d1 and d2 . . . . . 139LIST OF FIGURES1112131624313335535662657782889091946.3 Refractor elevation model constructed from accurate depth and dip samples 1416.4 Refractor elevation model constructed from inaccurate depth and dip samplesx1427.1 Subsurface first arrival wavefronts calculated with finite-differences7.2 Surface arrival time curves7.3 Finite-difference traveltime extrapolators7.4 First arrival wavefronts over a shallow syncine7.5 Reconstructed first arrival wavefronts over a shallow syndine7.6 Imaging a syndinal refractor7.7 Imaging an antidinal refractor7.8 Dependence of syncine on imaging time and downward continuation7.9 Reconstructed wavefronts and syncline image; noisy traveltimes7.10 Refractor velocity estimates for syncline and anticine models7.11 Shallow refraction traveltime data from Phalasarna, Crete .7.12 Shallow refractor image from Phalasarna data7.13 Comparision of observed and predicted traveltimes7.14 Forward and reverse first arrival wavefronts for final earth model148149151153154156157velocity . . 1581601611631641651668.1 Example raypath through cell ij of a 2D velocity model8.2 A 2D velocity model with wavefronts and raypaths from a surface source8.3 Five-cell constraint operator8.4 243 double-well VSP and crosswell raypaths8.5 Velocity tomograms obtained from VSP and crosswell arrivals; accurate data8.6 Velocity tomogram obtained from VSP and crosswell arrivals; noisy data8.7 81 crosswell raypaths8.8 Velocity tomograms obtained from crosswell arrivals onlyAl A noncentered, rotated effipseCl A simple 2D earth model with three layersGi Grid cells surrounding a source or receiver point176179• 183• 186• 187• . 188• . 190• . 191• . 204• . 2092191CHAPTER 1INTRODUCTIONThe traveltimes of head waves propagating in layered earth models have been studiedextensively since the inception of applied seismology in the 1920’s. The earliest Englishlanguage publications that examine this topic appear to be the classic works by Barton(1929) and Heiland (1929). These authors derive head wave traveltinie formulae appropriatefor 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 interfaceis horizontal, then head wave traveltime is a linear function of source-receiver offset distanceX:T(X) = mX + b. (1.1)The slope m and intercept b are independent of the recording geometry and are easilydetermined from the known earth model parameters. Numerous investigators, beginningwith Barton (1929), have extended this eicpression to multilayered one dimensional models.If the subsurface interface is clipping, then the same general mathematical formula applies. However, the slope and intercept are no longer invariant. Rather, the interceptdepends on the horizontal coordinates of the source (xs, ys) and the slope depends on theazimuth 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 directionof the refracting horizon. Hence, the azimuth angle EI is restricted to two values (say, = 0and ‘I’ = ir) that collectively refer to the updip and downdip directions. Surface-to-surfacehead wave traveltime isT(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 and2intercept 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 withplane, dipping interfaces by Ewing et al. (1939), Dooley (1952), Adachi (1954), Mota(1954), Ocola (1972), Johnson (1976), and Diebold and Stoffa (1981). The linear traveltimeexpression (1.2) also applies in this situation. Once again, the recording profile must beoriented 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 earthmodel parameters. Rather, each depends implicitly on the model parameters through a set ofraypath orientation angles. It is possible to determine these angles via repeated applicationof Snell’s law of refraction, starting at the critically refracting horizon and working upwardto the surface. Hence, equation (1.2) can be evaluated numerically.Methods for inverting observed head wave arrival times to obtain the earth model parameters are also discussed by the above authors. The techniques assume that traveltimesare recorded by inline forward and reverse spreads, or by a split spread. If the P-wave velocity of the uppermost layer is known, then the measured slope of a particular traveltimebranch can be used to infer the angle of incidence of the head wave at the surface. Repeatedapplication of Snell’s law of refraction from the surface downward determines the raypathorientation angle within the layer immediately above the critically refracting interface. Thedip and velocity of the critical refractor are then obtained from raypath orientation anglescalculated for forward and reverse propagating head waves. The method assumes that allvelocities and dip angles overlying the target horizon haye been previously determined. Finally, interface depths below the source locations are calculated from the measured intercepttimes.An alternative two dimensional inversion method exploits the variation of intercept timeb(x) with source position. Cunningham (1974) demonstrates that head wave arrival timesobserved on two inline spreads with the same source-receiver azimuth ‘P can be analyzed to3yield refractor dip, depth, and velocity. Although the slopes of the two traveltime curves arethe same, the intercept times associated with the different source locations provide sufficientindependent information for a solution. The method has practical utility in marine seismicexploration, where refracted arrivals are recorded on single-ended (i.e., unreversed) spreads.Cunningham (1974) examines simple models consisting of one and two layers overlying ahalfspace. However, it is possible to demonstrate the validity of the method for generalmultilayered models. The fact that both forward and reverse arrivals are not required for asuccessful 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 andinversion of head wave arrival times to three dimensional layered models. This class ofearth models is characterized by uniform velocity layers bounded by plane interfaces possessing arbitrary dip and strike. The linear traveltime expression remains valid in this threedimensional 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 throughoutthe interval 0 <2ir. Note that the slope in equation (1.3) depends only on the recordingazimuth pr• In general, the intercept depends on both the source location and the azimuthangle to the receiver. However, for the simple model consisting of a single layer overlying ahalfspace, 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 ofthe specified model parameters. Also, equations for critical and crossover distances arederived. All of these expressions are generalizations of more- familiar formulae that applyto 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 informationon the traveltimes of critically refracted waves in this basic three dimensional model. Dix4(1935), using geometric arguments, demonstrates that head wave isochrons (i.e., loci of fixedarrival time at the surface) are conic sections. A similar result is obtained by Dunkin andLevin (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 methodologyfrom 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-wavevelocities, two interface orientation angles, and a depth to the interface below a referencepoint. If the overburden velocity vi is known, then only four model parameters need to bedetermined 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 emanating 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 determine the attitude, velocity, and depth of the refractor. Two independent line profiles providefour measured data (two slopes and two intercepts) that allow an algebraic solution for thefour unknown model parameters. Inversion methods that are not based on conventionalline proffle recording geometries are also possible. For example, Chapter 2 indicates that allfive model parameters can be obtained by analysis of the geometric properties of head waveisochrons. This result suggests interesting possibilities for determining overburden velocityv from the refraction data alone.In Chapter 3, the derivation of the linear traveltime formula (1.3) is extended to multilayered earth models. A substantial simplification in the mathematical proof is obtained byusing a novel three dimensional form of Snell’s law of refraction. A previous proof given byDiebold (1987) is very cumbersome, and may actually be incorrect because it yields an expression that does not generalize to arbitrary source-receiver geometries properly. As in thetwo dimensional multilayered case, the slope and intercept in (1.3) depend on the orientationangles of the raypath within each layer. However, two angles are now necessary to describethe orientation of a raypath segment in three dimensional space. Chapter 3 describes a rapidcomputational procedure for obtaining these angles by applying Snell’s law repeatedly from5the 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 addressed in Chapter 4. Single layer and multiple layer algorithms are developed for dataacquisition geometries where sources and receivers are located on the surface. The inversion is posed as a constrained parameter optimization problem. An initial estimate of theearth model parameters is iteratively refined until an acceptable match is obtained betweenobserved and predicted arrival times. Similar parametric inversion schemes have been recently 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 andsurface-to--surface reflection traveltimes. Kanasewich and Chiu (1985) invert combined reflection and refraction arrivals. They use the iterative ray-bending approach of Chander(1977a) to calculate the traveltime derivatives needed for the inversion. The computationalprocedure developed in Chapter 3 is not iterative. Hence, it provides a rapid alternativeto ray-bending or ray-shooting methods for obtaining the necessary head wave traveltimederivatives.A novel feature of the inversion algorithms discussed in Chapter 4 is the introductionof constraint information via inequality relations that are imposed on the model parameters. Often, a priori geological or geophysical data are available to guide and constrain aninversion. For example, many shallow seismic refraction projects are undertaken in conjuction with a drilling program. Interface depth and layer velocity data may be obtained fromborehole logs. Constraints are particularly useful for the inversion of head wave traveltimes,because the problem can be very ill-conditioned and admit numerous solutions. Applicationof the constrained inversion algorithm to a crustal seismic refraction dataset from northernAlberta 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, thecentral portion of this thesis (Chapters 5, 6, and 7) develops specific improvements to varioustwo dimensional refraction traveltime inversion methods. By restricting consideration to two6dimensional models, examination of geological realities like undulating interfaces and/orvariable velocity media becomes mathematically tractable.Chapter 5 demonstrates that point estimates of interface depth and dip can be inferredfrom the observed refraction traveltimes. Two interpretation procedures are used for thediscussion. The Generalized Reciprocal Method (GRM) is a technique for delineating anundulating subsurface interface from refracted arrivals recorded on inline forward and reversespreads. It was developed by Palmer (1980, 1981) as an extension of the conventionalreciprocal method (Hawkins, 1961) of refraction traveltime interpretation. Although theGRM has been successfully applied to the problem of mapping undulating horizons, themathematical formulation of the method is based on a two dimensional earth model withplane interfaces. Hence, the improved head wave traveltime formula derived in Chapter 3can be applied to GRM analysis. A useful traveltime expression is obtained by specializingthe general equation to a two dimensional situation. As demonstrated in Chapter 5, all ofthe common GRM analysis tools can be derived from this novel 2D traveltime formula ina straightforward manner. Moreover, these new expressions allow point depth estimates ofthe refracting interface to be calculated from the measured traveltimes. A refractor depthprofile can then be obtained by interpolation. In contrast, the current GRM depth estimationprocedure 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 horizons into the model should yield more accurate results. Thus, Chapter 5 also proposesa head wave traveltime interpretation method tentatively named critical offset refractionprofiling. This inversion technique accommodates undulations in both the surface and therefracting 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 fromthe 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 dip7estimates obtained from critical offset refraction profiling provide additional constraints onthe interpolant.A detailed description of a method for calculating a smooth refractor depth profile froma set of point estimates of depth and dip is presented in Chapter 6. Determination of theinterface depth proffle is treated as an interpolation problem. Hence, the technique differssignificantly 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 forsolving this problem. The smoothest (i.e., minimum curvature) interpolant is the naturalmodel to adopt in this situation because refraction traveltime inversion methods assume,either explicitly or implicitly, that local interface curvature is negligible. Thus, the finalmodel for the refracting horizon is consistent with prior assumptions used for inferring itsdepth and/or dip. Posing the problem as an interpolation issue also has several specificadvantages regarding treatment of the data: (i) additional depth and dip values arisingfrom a variety of geological, geophysical, or engineering techniques are readily incorporatedinto 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 subsurfaceimage 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 pickedarrival times. The technique is an automated implementation of the classical WavefrontMethod for interpreting refraction traveltimes (Thornburgh, 1930; Hagedoorn, 1959; Rockwell, 1967). Modern finite-difference propagation algorithms are used to downward continuerecorded refraction arrival times through a near-surface heterogeneous velocity structure.Two such subsurface traveltime fields need to be reconstructed from the arrivals recordedon a forward and reverse geophone spread. The locus of a shallow refracting horizon is thendefined by a simple imaging condition involving the reciprocal time (the traveltime betweensource positions at either end of the spread). Refractor velocity is estimated in a subsequentstep by calculating the directional derivative of the reconstructed subsurface wavefronts8along the imaged interface. The principal limitation of the technique arises from impreciseknowledge of the overburden velocity distribution. This velocity information must be obtained from uphole times, direct arrivals, shallow refractions, and borehole data. Analysisof synthetic data examples indicates that the technique can accurately image both syncinaland antidinal structures. The method is also tested, apparently successfully, on a shallowrefraction dataset acquired at an archeological site in western Crete.Recently, tomographic techniques have been applied to the shallow seismic refractionproblem (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. Thefinal portion of this thesis (Chapter 8) presents an iterative tomographic inversion procedurefor reconstructing a two dimensional velocity field from measured first arrival times. Two keyfeatures of this technique are (i) use of a finite-difference algorithm for rapid and accurateforward modeling of traveltimes, and (ii) incorporation of constraint information into theinversion in order to restrict the nonuniqueness inherent in large scale, nonlinear tomographicinverse problems.Finite-difference traveltime computation (Vidale, 1988) provides a useful alternative toconventional raytracing in tomography. All first arrival wave types (direct and refractedarrivals, head waves, diffractions) are handled with relative ease. Curved raypaths, neededfor subsequent updating of the velocity model, are generated by following the steepest descentdirection through a computed traveltime field from each receiver back to the source. Sincearrival times are computed throughout a gridded two dimensional velocity field, very generalrecording geometries are easily accommodated. The main limitation of the method is that itis restricted to first arrivals. Hence, the algorithm developed in Chapter 8 cannot be appliedto 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, likesmoothness, on the constructed velocity model. In either case, the data equations are augmented with additional linear constraint equations and a least squares solution of the entire9system is obtained. In Chapter 8, the constrained inversion algorithm is applied to a simulated double-well VSP plus crosswell experiment. The results indicate that the algorithmcan accurately reconstruct a smoothly varying interwell velocity field. Inclusion of constraintinformation, in the form of horizontal and vertical first-difference regularization, allows thesolution of a traveltime tomography problem that is otherwise severely underdetermined(243 data and 10000 unknowns). For this example, no obvious improvement is obtained bythe addition of borehole velocity constraints.10CHAPTER 2HEAD WAVE TRAVELTIMES IN A THREE DIMENSIONALSINGLE-LAYER EARTH2.1 IntroductionMany refraction seismologists believe that a minimum of two reversed profiles, preferentiafly 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 belyzed to yield this same information. They state that similiar measurements made by onlytwo such profiles cannot define the three dimensional attitude of the dipping horizon. Thisstatement, however, is incorrect. Russell et al. (1982) do not fully utilize the informationcontained in the intercept times of the traveltime curves. One purpose of this chapter isto demonstrate that, in many cases, two refraction proffles are sufficient to define the threedimensional attitude, true velocity, and depth of a plane refractor. Generally speaking, themain condition required is that the two lines provide independent traveltime informationabout the subsurface.2.2 Earth model and recording geometryThe earth model consists of a single layer with P-wave speed v overlying a halfspace withP-wave speed V2. The plane interface separating the two media possesses, in general, a threedimensional dipping attitude. Let the zy plane of a right handed, rectangular coordinatesystem be coincident with the free surface of the layer; depth z is measured positive in thedownward direction. If r = xi + yj + zk is the position vector of an arbitrary point on thesubsurface interface, then the equation defining this dipping plane isrn = d. (2.1)11d is the perpendicular distance from the coordinate origin 0 to the plane and n is a unitvector pointing from 0 along the normal to the plane. Figure 2.1 illustrates that n isconveniently 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]ipdirection. 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).y0xnzFig. 2.1. Earth model and coordinate reference frame.12Figure 2.2a is a plan view of the surface recording geometry. The position vectors of thesource 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 ismeasured positive in the clockwise sense from geographic north. If the horizontal offsetbetween source and receiver is denoted by X (X 0), then the receiver coordinates can beexpressed 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 remainsfixed for a set of receivers that record energy from a single source.a) b)x xR0 0xSy yFig. 2.2. (a) Plan view of surface recording geometry. S and R denote a point source anda point receiver, respectively. (b) Spatial relationship of the two sources S1 and 52.132.3 Traveltime derivationConsider 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 raypathplane. In general, the raypath plane is not a vertical plane (parallel to the z axis). Only inthe particular case where the profile line is oriented directly updlip or downdip is the raypathplane vertical. The traveltime of a head wave propagating along the critical raypath fromsource to receiver can be worked out from simple geometric considerations. It isT =-+ (dS+dR) cosiV2 Vi(L Lcrit). (2.5)£ is the source-receiver range measured parallel to the refracting interface and d and dRare perpendicular distances from S and R to this interface. The critical refraction angle iis given by the usual expression: sin i = vi/v2. Of course, the head wave exists only forranges exceeding the critical range given byLcrit = (d + dR) tan i.R(2.6)Fig. 2.3. Head wave raypath. The plane of this diagram is perpendicular to the refractinginterface and thus is not necessarily vertical.S LdsdRP Q14The head wave traveltime formula is more useful to the geophysicist if it is expressedin terms of the horizontal offset X. A simple analytical derivation of the desired relation isgiven in this section. An advantage of this method is that it does not require visualizationof the three dimensional geometry of the problem. Starting with the general equation forthe dipping plane, it is easy to demonstrate that the perpendicular distances from S and Rto this plane are given by d = d — rs n and dR d — rj . n, respectively. HencedR = 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, theabove expression becomesdR = ds—X sinS. (2.7)It is stressed that the angle S is not the apparent dip of the refracting interface along theprofile 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 ismeasured in the raypath plane, whereas apparent dip is measured in a vertical plane. Forsmall values of true dip, S -y.The range £ is equal to the distance FQ in Figure 2.3. Hence15L2== (rS+dSn)—(rR+dRn) 2= (rS_rR)+(dS_dR)n2= Irs—rI?I2+2(dS—dR)(rS—rR)n+(ds—dR)== x2 — x2 s.ThusL = X cos8. (2.9)Equations (2.7) and (2.9) are the desired expressions. It is evident that these formulaecould have been derived by purely geometrical reasoning based on the raypath diagramof Figure 2.3. However, the relationship of the angle S to the proffle line azimuth and theorientation angles of the refracting plane would not have been easy to ascertain. Substitutingexpressions (2.7) and (2.9) into equation (2.5) and reducing yield an expression for the headwave 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)givesslnic /= 2d, . . 2.11cos(z—6)The crossover distance is defined as the particular source-receiver offset at which direct waveand head wave arrival times are identical. Setting T in (2.10) equal to X/vi and solving forX yieldscos zX035 = 2d . . (2.12)1— sin(z — 8)Fig. 2.4. Normalized apparent refractor velocity Va/V2. Each curve refers to a specific valueof the interface dip angle . The P-wave velocity ratio is vl/v2 = 3/5.16Finally, the traveltime formula (2.10) indicates that refracted arrivals propagate along thereceiver spread with an apparent velocity Va given byV1 Sifl 2Va . • V2. (2.13)sin(i — 8) sin(z — 8)The variation of apparent velocity with profile azimuth is depicted in Figure 2.4 for thespecific case where vl/v2 = 3/5.Equations (2.10) through (2.13) indicate that the three dimensional refraction formulaeare straightforward generalizations of those appropriate for the 2D problem: the angle Sreplaces the true clip S in the relevant expressions. It is emphasized again that S is not theapparent dip of the subsurface interface along the profile line. However, for a gently dippingrefractor, the practical difference between the two is small.3.02.5>>t51.00.50 100 200 300a—S (deg)172.4 Traveltime inversionTraveltimes recorded by a set of refraction proffles can be inverted to recover the 3Dattitude, true velocity, and depth of the refractor. Let the measured slope and intercepttime of the th such traveltime curve be m and r, respectively. From equation (2.10), theseare related to the assumed earth model parameters viasin(i — S) 2d. cosm = = , (2.14a,b)vi Viwith: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 linei. The perpendicular depth to the interface below source i is designatedHow many refraction lines are required in order to successfully invert for the earth modelparameters? Assuming that the overburden velocity vi is known (perhaps from boreholedata or traveltime analysis of the direct arrivals), then the earth model is defined by thefour parameters (v2, q, 9, d). Intuition suggests that two refraction profiles will yield the fourdata (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 8equations in 8 unknowns. Since the equations are nonlinear in the unknowns, a definitivestatement on the existence and uniqueness of a solution cannot be made. However, it canbe demonstrated by algebraic techniques that, in many cases, two refraction proffles aresufficient to solve the problem.The method proposed by Russell et al. (1982) requires three lines to obtain the refractor 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 casewhere the index i runs from 1 to 3, system (2.14) represents 12 equations in 10 unknowns;18hence, it would seem that redundant information exists. Although redundant data are always valuable in any practical inversion scheme, the theoretical minimum set of conditionsunder which the inversion is possible are of interest here.In the following analysis, only two refraction profiles are employed. First, define theangles 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 theclockwise sense from geographic north) and H is the horizontal distance between the twosource locations (see Figure 2.2b). Using these expressions, the angles 81, 8 and the depthsci, ds,, cl2 can be quickly eliminated from the system (2.14). A reduced system is obtainedconsisting 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 casesThree particular cases of the above system are examined before a more general solutiontechnique 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).19The two remaining equations contain three unknowns and cannot have a unique solution.Hence, two refraction proffles emanating from the same source location supply insufficientindependent information for a successful inversion. The split recording spread is a particularexample of coincident-source refraction profiles (Johnson, 1976).2) Parallel profiles. If the azimuths of both refraction lines are the same (cr1 = cr2), thenequations (2.17a) and (2.17b) above are not independent (identical azimuths imply identicalslopes of the traveltime curves and thus pi = pa). Hence, the parallel profile recordingconfiguration is not adequate to solve the problem either.3) Anti-parallel profiles. Suppose that line 2 is recorded in a direction opposite to that ofline 1 (cr2 = cr1 ± jr). Then equations (2.17a) and (2.17b) can be solved immediately for thecritical refraction angle:= 1U1 + P2 (2.18)The refractor velocity is obtained via V2 = v1/ sin c. This procedure for calculating V2 isidentical 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) recordedover a three dimensional dipping interface.With the critical angle i determined, the system (2.17) can be solved for the remainingunknowns. The azimuth of line 1 is designated by a. Straightforward, but tedious, algebraicmanipulation then yieldstanO —— (r—7-2) cosa + H(mi —m2) cos/3 (219)—(Ti — T2) sin a + H(mi — m2) sin/3sin 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 byF..(2.21)2H cos i sin(a — /3)20A minor ambiguity associated with selecting the proper branch of the inverse tangent function in equation (2.19) is easily resolved by ensuring that the angle 8 satisfies the originalsystem (2.17). After 0 and çb are determined, one of the equation pairs (2.14b,d) can besolved 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) containsan important indeterminacy when the azimuth angle /3 of source 2 equals a or a ± ir. In thiscase, the recording geometry consists of two colinear, anti-parallel proffies. The identity ofreciprocal times recorded at the two inline shot positions (i.e., mlH+Ti m2H-Hr)impliesthat each equation reduces to 0/0. Hence, two colinear anti-parallel proffles are inadequatefor defining three dimensional refractor attitude. The reversed spread is a typical exampleof this configuration. As long as there exists some perpendicular offset distance betweenthe two lines recorded in opposite directions, the inverse problem is well posed, at least in atheoretical sense. The latter data acquisition geometry is a fairly common arrangement formany land and marine seismic surveys.Arbitrary line azimuthsA 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 treatedin the previous section are excluded. Eliminating 8 and çS results in a quadratic equation intan i:A tan2i+ B tani + C = 0. (2.22)The coefficients A, B, and C depend upon measured quantities and recording geometry asfollows:A V1T1T2 f= 2Hsina —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)21In the case of anti-parallel proffles the constant A vanishes and tan i is obtained bysolving a linear equation. It is easy to verify that the expression for the critical angle igiven 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)2Awhere the sign of the radical has been chosen to yield an indeterminacy for A = 0. Application of L’Hopital’s rule then indicates that (2.18) is recovered as A —* 0. After iis determined, simply back substitute sequentially into the expressions for the remainingunknowns. This is done numerically rather than symbolically since the formulae rapidlybecome unwieldly. Hencetan9 =— H cosi sin(ic — p.i) cos/3 — [vl(T1 — r2)/2] C05a1 (225)H cosi sin(i— i”) sin3 — [vl(r1 — ‘r2)/2] sinaiAn 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 theoriginal system of equations. Any one of the expressions (2.17a,b,c) can then be solved forthe true dip angle q5. Finally, depth d is determined as previously described.2.5 Traveltime isochronsHeretofore, the analysis is concerned with an inversion method for critically refractedtraveltime data observed on two line profiles. An alternate approach to interpreting threedimensional refraction data entails recording arrival times on an areal grid of receivers. Themeasured data are posted and contoured on a map display of the receiver array. Characteristic patterns observed in the contoured data are used to qualitatively infer subsurfacestructural style. If the traveltime data are sufficiently precise, a quantitative estimate ofvarious structural or stratigraphic parameters (dip and strike angles, locations of faults or22edges, magnitude and sense of fault throw, refractor velocity, etc.) is possible (Dunkin andLevin, 1971; Palmer, 1986, p. 246-249).For a fixed traveltime T, the head wave wavefront occupies a three dimensional spatiallocus consisting of a frustum of a cone. The intersection of this locus with the recordingplane z = 0 is defined to be a traveltime isochron. A closed form expression for a headwave 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 thedistance 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 Appendix A. Hence, a locus of fixed head wave arrival time T is a conic section (either anellipse or a hyperbola) with principal axes parallel to the dip and strike directions of therefracting interface. Furthermore, the parameters defining the geometry of the curve canbe 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 qSe — , a —. 2 .2 — • 2 .2sin z sin z — sin sin z — sin(2.28a, b, c)These three expressions are consistent with equivalent formulae derived in an entirely different manner by Dix (1935). Two distinct cases must now be examined.23Case 1: > q1. If the critical angle exceeds the interface dip angle, then the eccentricityof the conic section is less than one (e < 1). Hence, traveltime isochrons are effipses thatencompass the source position. However, the shotpoint is not located at a center or a focusof any effipse. Rather, the center is displaced in the updip direction by the distance I, andthe 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) withradii that increase linearly with T — T.Case 2: i < . In this situation, the eccentricity exceeds one (e > 1) and isochrons arehyperbolae. This case may arise where there is a strong velocity contrast between overburdenand refractor, together with a large dip on the intervening interface. Head wave traveltimesactually decrease with offset distance in a wedge of azimuths centered about the updipdirection (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-receiveroffset 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 withequation (A4), the parameters defining this conic aresin 4’ 2d sin i, cos cos 4’ —2d5 sin2 i sin 4’e = , a= 2 • 2,= 2 • 2cos i cos z — sin 4’ cos i — sin 4’24400200E-cc—20—4004002000-c0c —200—400—400 400eastingFig. 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, andthe crossover distance by a dashed curve. The source (asterisk) is located at the origin, andthe straight line indicates the outcrop locus of the refractor. (a) Earth model: vi = 1500m/s, v2 3000 m/s, q. = 15°, 0 = 45°, l = 100 m. Isochron eccentricity = 0.518. (b) Earthmodel: v = 1500 m/s, v2 = 5796 m/s, = 30°, 0 = 45°, h = 100 m. Isochron eccentricity= 1.932.0(m)25400200E 0-c-—200—400400(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 ellipsesurrounding 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, lowvelocity contrast environment where rays critically refracted downdip never return to thesurface.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 realtraveltimes. Figure 2.5a illustrates elliptic traveltime isochrons ( > 4) and an ellipticcritical 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 0easting26elliptic (ia> q) and the critical distance curve is hyperbolic (i + q> Tr/2). The remainingcase (hyperbolic isochrons and hyperbolic critical offset) is not illustrated. The crossoverdistanceX03(a), computed from equation (2.12), is depicted by a dashed curve in eachpanel. Since expression (2.12) cannot be put into the form of equation (A4), this curve isnot a conic section.2.6 Inversion of traveltime isochronsIf the geometric parameters of a traveltime isochron can be determined from head wavearrival times observed on an areal recording array, then it is possible to invert for the earthmodel parameters. The principle axes of an effiptic or hyperbolic isocliron are parallel to thedip 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 termsof e, a, and yield1 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; theargument of each square root is still positive.Equations (2.30a,b) indicate that the critical and dip angles can be determined fromthe geometric parameters of an observed traveltime isochron. If the overburden velocityvl is known from ancillary data, then(2.30c) can be solved for the intercept time r. Thenormal 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 respectto recording profile azimuth c. Thus, a value for T can be obtained by linearly extrapolating27to 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 v1gives__________________—e2)(a —l)Vj= eIT—rI.(2.31)Hence, the overburden velocity can be calculated from the refraction data alone! With v1determined, the refractor velocity is obtained via v2 = vi! sin i. Interestingly, the abovesolution 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 indeterminate. Overburden velocity v must then be obtained from independent measurements.Practical issues related to fitting a conic section to a finite set of spatially sampledarrival times are not addressed here in any detail. Rather, the intention is to demonstratethe 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 contaminatedtraveltimes 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 maybe necessary. Finally, it may be possible to fit several conics with different isochron timesT to the measured data. Statistically robust estimates of the earth model parameters canthen be obtained by using some form of averaging.2.7 ConclusionAlthough the earth model examined in this chapter is quite simple, there are undoubtedly 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 propagationwithin the model yields closed form mathematical expressions for traveltime, critical andcrossover distances, intercept time, and apparent refractor velocity. These expressions aregeneralizations of more familiar formulae that apply in one and two dimensional situations,and are valuable for forward modeling purposes.28Algebraic solutions of the traveltime inverse problem for the same earth model are alsoderived. There are numerous practical recording geometries where the analysis of criticallyrefracted arrival times recorded on only two unreversed spreads can yield the three dimensionai attitude, true velocity, and depth of a plane subsurface horizon. The well knownmethod for determining refractor velocity in a 2D situation (measurement of apparent velocities on colinear forward and reverse profiles) applies directly to the more complicated3D model. However, it is not necessary that the proffles be colinear; anti-parallel profilesare 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 multilayered earth models will be valuable for defining the minimum requirements for a successfulinversion of three dimensional refraction traveltime data.Finally, study of head wave arrivals recorded in an unconventional data acquisition geometry (i.e., an areal receiver array) reveals an interesting possibility for estimating overburdenvelocity from the refraction data alone. At present, this is primarily a theoretical result, butit does provide some justification for the numerical inversion results obtained in Chapter 4.29CHAPTER 3HEAD WAVE TRAVELTIMES IN A THREE DIMENSIONALMULTILAYERED EARTH3.1 IntroductionThere is a notable paucity of papers on the subject of head wave propagation in threedimensional, multilayered earth models. Chander (1977b) examines a model consisting ofuniform velocity layers separated by plane interfaces with arbitrary strike and dip, and describes a method for calculating head wave traveltimes between specified source and receiverpositions on a horizontal surface. If an array of receivers is colinear with the source, the headwave 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 straightline. Chander (1977b) locates the two initial points via raytracing techniques.Chander’s work is purely numerical and does not provide much insight into the dependence of head wave traveltime on the parameters that define the earth model. Moreover, itis restricted to conventional data acquisition geometries. Buried sources and/or receivers aswell 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 threedimensional multilayered earth, but derives traveltime formulae for reflected and criticallyrefracted waves. These formulae are logical extensions of familiar traveltime expressionsthat are appropriate for one and two dimensional layered models. Thus, they offer thepossibility for extending several known traveltime inversion techniques to accommodate 3Dplanar structure. Unfortunately, Diebold’s derivations are very ambiguous. Furthermore, hisgeneralization 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 becredited with an original contribution to traveltime analysis for this particular class of earthmodels.30This chapter provides a rigorous derivation of the three dimensional head wave traveltimeformula. A related expression for the traveltime of a reflected wave propagating in the sameearth model is obtained as a byproduct of the analysis. The mathematical proofs of theformulae 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 arbitrary3D recording geometries and/or mode-converted waves. Finally, a rapid numerical methodfor computing the arrival times of critical refractions is presented, and is illustrated withsimulated examples from shallow refraction exploration and vertical seismic profiling (VSP).3.2 Earth modelConsider an earth model consisting of a set of homogeneous and isotropic layers boundedby plane interfaces. In general, each interface may possess a three dimensional dippingattitude. 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 inthe downward direction. The locus of plane interface i satisfies the equationr•n = dj, (3.1)where n is a unit vector normal to the interface and dj is the perpendicular distance from0 to the interface. Figure 3.1 indicates that n is conveniently described by two interfaceorientation 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 coordinates31y0Fig. 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 ofits Cartesian components:n = nj,i + ni,,i +xnizvi +1with n + m, + = 1. This convention is followed in the sequel.32Solving equation (3.1) for z as a function of z and y yields the vertical depth of interfacezQc, y) z(O, 0)——y [3], (3.2)Th:,zwhere 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 arenumbered sequentially in the downward direction. Interface i overlies layer i. The verticalthickness of layer i is defined to be h(z,y) zj+1(x,y) — zj(x,y). Thush(x,y) = h(0,O) + [Th—hh+ij— (3.3)‘i,z i+1,z fl,z i+1,zwhere h(0, 0) z+i(O, 0) — z(O, 0) is the vertical thickness of the layer beneath thecoordinate origin 0.Finally, the seismic wave propagation speed assigned to layer i is given by v. This maybe either the compressional wave speed cx or the shear wave speed f3. This fiexibffity allowsthe resulting traveltime equations to apply either to P, S, or mode-converted waves.3.3 Raypath geometryInitially, the analysis is restricted to the case where both the source and receiver arelocated on the surface. Generalization to an arbitrary data acquisition geometry is straightforward and is given in a later section. The horizontal coordinates of the point source S andpoint receiver R are (as, Ys) and (R vi?), respectively. Their vertical coordinates are easilyobtained 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 isdivided 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. Thepropagation time along each portion is calculated, and then all three are summed to obtainthe surface-to-surface head wave traveltime.33Fig. 3.2. Schematic representation of the raypath of a wave critically refracted on interfacek 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 portionof the raypath, the propagation direction within layer i of the wave critically refracted atsubsurface interface k is described by the unit vector pj:Pik = Pik,z + PiJc,.i +SR....P0Similarly, the upward propagation direction within layer i of the wave critically refracted34from 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, = qkkAt 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 depictedin Figure 3.3a. The plane of this diagram is the plane of incidence defined by the incidentpropagation direction Pi—1,k and the interface normal n. Snell’s law of refraction consistsof the following two conditions:(i) the unit propagation vector of the transmitted ray (pjk) is contained in the plane ofincidence,(ii) sin Pjvj_1 = sin v/v, where and v are positive acute angles measured from the interface normal to the incident and transmitted propagation directions Pi—1,k and Pik,respectively.Both conditions are contained in the single vector equationfli )< Pi—i,k = n X Pik (3 4)Vj_1 VThe vector formed by these cross products points out of the plane of the diagram in Figure3.3a. In component form, equation (3.4) is35Fig. 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-1Pi-1,k aInterface IniP1kiflterfa Iqi-1,k bv’—1qikni361 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_ Vtii 1 ir 1[fli,zPi—1,k, — fli,vPi—1,k,zj — [fli,zPik,i, — flz,Pik,zj. (3.5c)Vj_1Similarly, when the upgoing wave encounters interface i from below, Snell’s law in theformflj X q—1,k ni x (3 6)Vj_1 Vholds (see Figure 3.3b). In this case, the angles 1ti and v exceed ir/2 radians. The componentform of expression (3.6) is analogous to equations (3.5).The three dimensional statement of Snell’s law of refraction given by the above expressions 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 expressionsare equivalent to the raytracing formulae (equations (3.27) and (3.28) below). The value ofthe current formulation is that it leads to a substantial simplification in the mathematicalproof of the traveltime equations.3.4 Traveltime derivation3.4. 1 Downgoing traveltimeAn expression for the traveltime increment of the downgoing wave as it traverses layeri is derived first. Let the position vectors r: and r:+I denote the intersection points of thedowngoing ray with interfaces i and i + 1, respectively. Then rj+1 = r + 13 Pik, where374 is the length of the (straight line) raypath segment within layer i. Solving for 4 gives4 Pk (r+I — rj. In component form, this expression is4 = 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 layervelocity v:ii Pik,z Pik,y Pik,zti = — = ——(xi+i — x) + —(yi+1 — yi) + —(zj+i — z).vi v vi viThe vertical (z) coordinates of the intersection points can be expressed in terms of thehorizontal (x, y) coordinates by using the equation for a dipping plane interface. Equation(3.2) yieldsrni,zl F’i,ylz zj(xj,yj) = zj(O,O)—aji—-—i —yi——--,LThi,zJ Lfli,zJ—rni+1,zl rni+i,yZj+1 = zj+1(xj+1,yj+l) = zi+1(OO)—xi+1[ ] Yi+4Hence1i,1 f2i,v1 1’i+l,z1 1’i+i,y1—= h(O, 0) + z—j + — x+1 [ j — Yi+1 [Th:,zwhere h(0, 0) = z+i(0, 0)— z(0, 0) has been used. Substituting this result into the equationfor 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,zvi -I L v38The total downgoing traveltime is obtained by summing all of the individual traveltimeincrements 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 VVThe sum involving aji and Yi-i-1 is now reindexed and combined with the other sum. Theresult isk—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=2Vj_1 vi39At 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! Thedowngoing traveltime from (x1,yl) = (xs,ys) to (zk,yk) (xp,yp) reduces toIc—i___________.Ldownvi+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)Z3.4.2 Upgoing traveltimeThe traveltime along the upward propagating portion of the total raypath is derived bysimilar techniques. Snell’s law in the form (3.6) is used at each interface in the overburden.The result isk—i‘ç—___‘up— Ldi=ivi—XQ(flk,zqk.i,k,z — flk,zqk_1,k,z) + yQ(nk,zqk_i,k,y — nk,yqk1,k,zVk.lflk,z+ZR(Th1,zqik, — ni,zqlk,z) + yR(ni,zqlk,y — n1,qik,z) (38)v’n’,z403..3 Critically refracted traveltimeThe final traveltime increment needed for the derivation corresponds to the criticallyrefracted segment of the total raypath. Let the position vectors rp and rq refer to theintersection 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 twopoints (and thus lying entirely within the plane of interface k). Then rq = rp + tk Pick andhence l = Pick . (r — rp). The propagation time along this path length isPlck,z Pkk,y Pkk,ztic = —=(x — zp) + (y—yp) + (z — zp).Since points P and Q reside on the same plane interface, equation (3.2) yieldsrk,z1ZQZp = Zk(Q,yQ)-Zk(Zp,yp) = —(q---a)-——jfk,z k,zThe 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 traveltimeThe total surface-to-surface traveltime of the wave critically refracted on interface k isobtained 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) yields41k—iTk(XS,YS,aR,YR)— qiie,)vii=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 refractionand is given byF(xp,yp,XQ,yQ) =xpF 1 1I (k,zPk—1,k,z — Thk,zpk—1,k,z) — (‘k zPkk,z — ThkzPkkz)]k,z Lvk_1 VkFi 1+ I (flk,zPk—1,k,y — flk,yPk—1,k,z) — (Thk,zPkk,y — flkPkkz)]k,z [vk_1 Vk_1I (n,zq_1,,z — nk,zqk_1,k,z) — — ThkzPkkz)]k,z [Vk_i Vk1I (‘k,zqk—1,k, — nk,yqk_1,k,z) — (k,zPkk,y — flkiPkkz)].k,z [vk_ISince the wave is critically refracted at interface k, the propagation direction vectors Pkkand kk are identicai. Then, requiring Snell’s law to be satisfied at points P and Q results42in F(xp, yp, XQ, y) = 0. The final formula for the surface-to-surface traveltime of the headwave becomesk—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 formulaObviously, if the source or receiver is located at the coordinate origin, then the traveltimeformula (3.10) simplifies considerably. Another simplification arises with a. horizontal surface(ni, = = 0, n = 1). Expression (3.10) reduces toTk(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 specifyingthe 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 areXR x+Xcos’I’, YR = ys+Xsin.43Substituting these expressions into equation (3.10) yieldsk—iI V T ç’ h(0, 0)(Pk,z — qik,z)1kkXs,yS,,) =Thi,z(pik,x— qik,z) — 21,z(P1k,z — qlk,z)— 2sViThi,zfll,z(pik,y — qlk1y) — Th1,y(Pik,z — qik,z)—ysVifl.1,z+cos (n1,zq1, — Th1,zqlk,z) + Sfl (ii,zqiic,y — ni,yqik,z) (3.12)Vlfli,zNote that X is the horizontal distance between source and receiver; the actual distance islarger since it is measured within the plane of interface 1. It is straightforward to demonstratethat the true source-receiver distance is L XsJl + tan2 li cos2(— 8), where 4’i and 6are the dip and azimuth angles of the surface.Equation (3i2) is an extension of the common ‘slope and intercept’ traveltime formulafor head waves to three dimensional, multilayered earth models. For the particular case ofa model with only two layers and a horizontal surface, it can be shown that this equationreduces to the closed form expression (2.10) derived in the previous chapter. This serves asan important check on the correctness of the general formula. However, the proof entailssome cumbersome algebra, and hence is relegated to Appendix B.Finally, consider the specialization of the general traveltime formula to a two dimensionalearth model. In the model parameterization utilized here, all of the y components of theinterface normal vectors vanish (equivalently, all interface azimuth angles are restricted to= 0 or = ir). Additionally, the recording profile must be oriented perpendicular to thestrike directions of the subsurface horizons. Hence, I1 = 0 or r also. If this second condition44is not satisfied, the wave propagation vectors p,j and qjj are not confined to the xz planeand a full three dimensional treatment is necessary.Since many previous investigators have assumed a horizontal surface, equation (3.11) isused as the point of departure for the analysis. Setting s = YR = 0 givesTk(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 equationsdeveloped by Diebold and Stoffa (1981) and Diebold (1987). If the source is located at thecoordinate origin, then (3.13) is also consistent with earlier two dimensional head wavetraveltime formulae published by Dooley (1952), Adachi (1954), Ocola (1972), and Johnson(1976). All of these investigators prescribe vertical layer thicknesses, either beneath theorigin or the shotpoint. In contrast, Ewing et al. (1939) and Mota (1954) measure thicknessnormal to the basal interface bounding a layer. Hence, their traveltime equations, althoughdesigned to treat an equivalent situation, differ in mathematical detail.3.5 Generalizations of the traveltime formulaHeretofore, both the source and the receiver have been restricted to the surface. Moreversatile formulae are needed to model data acquisition geometries with buried sourcesand/or receivers. These situations arise in surface-to-borehole, borehole-to-surface, andborehole-to-borehole seismic experiments, as well as with placement of sources and/or receivers in underground mines.3.5.1 Source and receiver on separate interfacesLet the source S be located on the th interface with 1 < k. The downgoingtraveltime is obtained by summing the layer traveltime increments t, i = j,j + 1,. . . , k — 1.45Equation (3.7) generalizes tok—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,zSimilarly, if the receiver R is located on the 1 interface (1 1 < k), then the upgoingtraveltime 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,zThe critically refracted traveltime increment is still given by equation (3.9). Summing thisexpression together with (3.14) and (3.15) yields the total traveltime:46h(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,zwhere zs = zj(xS,yS) and ZR = zl(rR,yR). This is the proper expression for head wavetraveltime when source and receiver are located on different interfaces of the model. Itdiffers significantly from the analogous formula published by Diebold (1987). His expressionis actually a special case of the general equation (3.16); in particular, it is only valid if both thesource 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 examinationof a simple two dimensional earth model in Appendix C. The analysis demonstrates thatequation (3.16) reduces to the known traveltime solution for this situation, whereas Diebold’sexpression yields an erroneous result.3.5.2 Arbitrary source and receiver locationsA further generalization is obtained by allowing the source and receiver to be locatedwithin designated layers. Assume that the source is located in layer j at a vertical depthds below the the immediately overlying interface (the ih). Similarly, let the receiver belocated within layer 1 at a depth dR beneath interface 1. These incremental source andreceiver 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 thissituation. Traveltime increments induced by the source layer j on the downward path andthe receiver layer 1 on the upward path must be treated separately, because the wave doesnot propagate across the full thickness of each layer. The result of the analysis is47h(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,zwherezs = zj(xs,ys)+ds and ZR = zz(XR,yR)+dR. Note that the prior expression (3.16) isrecovered in the limit as — 0 and dR —* 0, as expected. As an additional check, examinethe 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 berewritten in terms of the thicknesses at the source and receiver positions (h(xs, ys) andhl(XR,yR)) via expression (3.3). Then, by appealing to Snell’s law at each basal interface,equation (3.17) is recast ashi(0, O)Pilc,z [ha(xs, ys) — dsjpk,Tk(XS,YS,ZS,XR,YR,ZR)=.+1=2+17)3—h(O,0)q,—[hz(XR,yR) —vi vjizzj+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,z48Comparing equations (3.18) and (3.16) indicates that the proper form of the traveltirneexpression 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 thetraveltime formulae: asymmetric wave propagation paths can be treated. An asymmetricraypath is defined as one where the mode of dowr&going wave propagation in the ith layerdiffers from the mode of upgoing wave propagation across the same layer. Strictly, different symbols should be used to designate the wave speeds within layer i in the downwardand 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 maintain notational simplicity. The velocity v appearing in each sum is simply interpreted asthe 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 wavespropagating in a three dimensional layered earth model.3.5.3 Arbitrary reference points for layer thicknessIndividual layer thicknesses enter the traveltime expressions evaluated at the coordinateorigin 0. An alternate form of the traveltime equation is characterized by layer thicknessesspecified below the source and the receiver. This variant is particularly suitable for the timeterm, delay time, and reciprocal time inversion methods. Hence, the previous derivation isnow modified to incorporate layer thicknesses prescribed at arbitrary reference locations;these points can then be specialized to the source and receiver positions. The resultingtraveltime expression forms the point of departure for a three dimensional extension of theaforementioned inversion techniques.The depth of the i interface, referred to an arbitrary location A with horizontal coordinates (XA, VA), is given byri,y1zj(x,y) = zj(XA,yA)—(x—xA)[—--—j—(ii— YA)[zj. (3.19)n:,z ni,z49Consider the downgoing raypath first. The previous expression for the traveltime increment Lj induced by wave propagation across the th layer is trivially modified toPilc,z r 1 Pik,y 1 1 Pik,zt =____—ZA)—(x XA)j +— YA) — (Yi — YA)j + ——(z+ — z).vi vi viThe vertical (z) coordinates of the ray intersection points are now expressed in terms of thecorresponding horizontal coordinates by using equation (3.19):IThizi Ifliyzj — z = h(xA,yA) + (x— X4)Lj + (y —lil,z n:,zlmi+1,zl______—(x+ — XA)—(y+i — YA) I I’Lflj+1,zJ Lfli+1,zJwhere h(xA,yA) zj+1(XA,yA) — zzfrA,YA) is the vertical thickness of layer i below thereference 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,zVi—(x—XA) i,zPik,z — i,zPik,z — (Vi — VA) i,zPiIc,y i,yPik,zni,z vi ni,z viThe traveltime increments t are now summed over all layers in the overburden. ApplyingSnell’s law at each plane interface results ink—is—.,___________1down =Vii:=i50+(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 ,zA similiar analysis yields the upgoing traveltime. Layer thicknesses are now referred toa 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,zFinally, the critically refracted traveltime increment is given by a simple alteration to theprior 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,z51The total head wave traveltime is now obtained in the usual manner by summing thecontributions along all three raypath portions. Adding equations (3.20), (3.21), and (3.22)givesk—ih(zA,yA)pk, — h:(XB,YB)q:k,zTk(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 atthe critically refracting horizon demonstrates that F vanishes identically. Furthermore, itis common (but not mandatory) practice to select the reference points A and B to becoincident 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 dramaticallytok—ih(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,2The two dimensional form of equation (3.23) has some similarity to an analogous expression published by Ocola (1972), in that individual layer thicknesses are measured vertically52below the source and receiver. However, it is quite different from Palmer’s (1980) 2D trayeltime equation, which currently forms the theoretical basis of the generalized reciprocalmethod (GRM) of seismic refraction interpretation. Following Ewing et al. (1939), Palmermeasures thickness perpendicular to the base of each layer, and thus arrives at a differentexpression. Nevertheless, the GRM can be formulated on the basis of the two dimensionalversion of equation (3.23); a distinct advantage is then gained in constructing the refractordepth proffle (see Chapters 5 and 6).- 3.5.4 Reflection traveltimeThe previous analysis has concentrated on. waves that are critically refracted at thekth subsurface horizon. However, it is straightforward to demonstrate that the traveltimeformulae also apply to waves that are reflected from interface h of the model. Snell’s lawof reflection for this situation is illustrated in Figure 3.4 and is compactly expressed ask X Pk—1,k — k X qk—I,kd — uVklNote that equation (3.24) allows for a possible mode conversion upon reflection, i.e. v1 isnot necessarily equal to v_1.For a wave reflected at interface k, points P and Q are coincident; there is no interveningcritically refracted raypath segment. The total traveltime is obtained by simply adding thedowngoing and components. Hence, setting (XQ, yq) = (xp, yp) and summing equations(3.14) and (3.15) yields53Fig. 3.4. Reflection of a raypath at interface k of the earth model. If a mode conversionoccurs, 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 pointPk-1,k qk-1,kkand is given by54G(xp,yp) =xp 1 1d (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—iyp 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 13klDistinct downgoing and upgoing layer velocities are explicitly incorporated into the aboveexpressions 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 derivedfor head waves!3.6 Rapid traveltime computationIn order to compute traveltimes via the above formulae, the unit propagation vectors Pikand qjk overlying the refracting/reflecting interface must be determined. For the reflectionproblem, this set of vectors depends on both the offset distance X and the azimuth angle 1J1of the receiver relative to the source. However, the critical refraction problem is qualitativelydifferent; the propagation vectors depend only on the azimuth ‘I’. This particular featurecan 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 asPk(’P) and qik(), although the explicit dependence on ‘P is often suppressed for notationalconvenience. The functional form of this dependence is not known. However, with a minimalamount of raytracing, it is possible to numerically generate the function ‘P(pjk, q) over thefull range of possible recording azimuths (2z radians). Inversion of this function then yieldsthe propagation vectors for a prescribed value of the source-receiver azimuth angle. Thistechnique is discussed in both general and mathematical terms in the following two sections.553.6.1 General descriptionThe computational algorithm is based on the close relationship between a criticallyreflected and critically refracted raypath. For a given source-receiver azimuth angle, bothraypath 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 reflected/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 orientationof this vector within the plane of interface k is defined by an angle x measured from anarbitrary 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) relativeto the interface normal and is contained in the plane defined by Pkk and Theappropriate three dimensional form of Snell’s law is applied at all interfaces interveningbetween 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. Theincidence angle of the arriving vector pk1,k (in the plane defined by Pkk and nk) alsoequals the critical angle. A three dimensional ‘backward propagating’ form of Snell’slaw 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 xhas been incremented by a total of ir radians.56.. ..Fig. 3.5. Schematic representation of a raypath critically reflected/refracted at point P oninterface k. S and R denote a surface source and receiver, respectively.Pk--i,kqk-1,k57This procedure numerically defines a function = f(x) over an interval {xo, xo + it]. It isnot 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. Thissymmetry condition arises from raypath reciprocity: reversing the direction of the criticallyrefracted propagation vector Pkk merely interchanges the positions of the source and receiveron the surface.If the propagation vector Pkk makes one complete rotation on the critically refractinginterface (i.e., x increments by 2ir radians), points R and S make one complete closed circuiton 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 appropriatevalue of x for a specified source-receiver azimuth angle. Finally, unit propagation vectorsPik 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 azimuthscontained in the data acquisition geometry are treated by this same function. The numericalmethod for inverting g(x) is briefly stated in the next section.3.6.2 Raytracing techniqueThe pertinent mathematical details of the aforementioned computation procedure arenow described. A distinction must be made between the downgoing and upgoing portionsof the critically reflected raypath. Thus, superscripts ‘d’ and ‘u’ refer to these raypathsegments, respectively.The Cartesian components of the critically refracted propagation vector Pkk are chosento bePkk,z = cosXcoscbkcosOk — sinxsin, (3.26a)= cos x cos sin 0k + sin x cos 6k, (3.26b)Pkk,z = — cossinq. (3.26c)58It is straightforward to verify that Pkk 1 and Pkk = 0. Hence, Pkk is a unit vectorperpendicular to the interface normal k• Its orientation within the plane of interface k isdetermined by the angle x If x 0, then Pkk points in the maximum updip direction alongthe critically refracting horizon. Increasing the angle x rotates Pkk in a clockwise sense whenviewed from the coordinate origin 0.Consider the downgoing portion of the raypath first (segment SF). Shah (1973) usesthe following three dimensional form of Snell’s law for the ray refracted at interface i:d V1 dPa —P—i,k + (cos ii — — cos i )n1.Vj_1 V_Hence, the transmitted propagation vector is a linear combination of the incident propagation vector and the interface normal vector. Solving for Pi—i,k yieldsPi—1,k = + (cos ,.d — cos vd)n1. (3.27)Vt VtThis form gives the incident vector in terms of the transmitted vector and the interfacenormal. The cosines of the incidence angles can be obtained by the following four-pointprescription: (i) cos = Pa . n, (ii) sin d v”l — cos2 vd, (iii) sin = (v_1/ )sin .A,(iv) cos /1 — sin2Similarly, the upgoing ray (segment FR) at interface i satisfies Snell’s law in the formii U Uqi—1,k = —--qik + (cosv ————cos/1 )n, (3.28)Vt Vtwhere 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 anegative sign is used to calculate cos v since the angle v exceeds 7r/2 radians. Startingwith pj (= kk) given by expression (3.26), equations (3.27) and (3.28) are evaluatedrecursively for i = k, k — 1, k — 2,. . . , 3, 2 to generate the unit propagation direction vectorsin all layers overlying the critically refracting interface.59Intersection points of the critically reflected raypath with the plane interfaces of themodel are easily determined after the propagation directions are known. For the downwardpath segment within the i’ layer, these are related viad dr2= —3.29Hence, 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)) yields1d — (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 evaluatedrecursively with i = k — 1, k — 2,. . . , 2, 1 to yield the ray intersection points r and r onthe surface. Finally, the azimuth angle E’ of r (the receiver point) with respect to r (thesource point) is calculated by= tan1 — . (3.31)—iJ is obviously a function of the reference angle x: ‘I’ = g(). In principle, the inversefunction x = g1() can be determined. In practice, numerical calculation of the complete inverse function g1 is unecessary. Rather, inverse interpolation between neighboring60computed values of ‘I’ is used to locate the particular value corresponding to a specifiedsource-receiver azimuth angle ‘P. Linear interpolation is adequate, provided that the sampling in x is sufficiently fine. The value is then used to regenerate the complete set ofpropagation vectors Pik and q needed to evaluate the head wave traveltirne formula.A comparison with an example presented by Chander (1977b) validates the accuracy ofthe calculation. He considers a four layer model defined by the parametersci = 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 criticallyrefracted on interface 4 are determined for a recording azimuth ‘P = 33.69° (the samplinginterval used for x is 0.25°). In particular, the downgoing propagation vector at the source‘SP14(33.69) = (0.523126)i + (0.322360)j + (0.788938)k.Chander (1977b) iteratively solves two coupled nonlinear equations for the raypath takeoffangles at the source, given the above azimuth angle to the receiver. The orientation of thetakeoff vector P14 can be described by two angles in a manner analogous to the interfacenormals (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.613.7 Forward modeling examplesTwo examples presented in this section illustrate the utility, as well as some of thelimitations, of the head wave traveltime formulae for forward modeling applications. Therecording geometry in the first example is the common reversed proffle; all sources andreceivers are located on the surface. The second example considers a typical offset VSPgeometry (surface source and downhole receivers).3.7.1 Profile geometryEquation (3.12) expresses the head wave traveltime in terms of the horizontal offsetdistance between a surface source and a surface receiver. It is written in condensed form asTk(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 expression is a three dimensional extension of the slope/intercept formulae that are commonlyused to describe head wave traveltimes. Equation (3.32) is evaluated for a shallow three-layermodel defined by the parametersqSi 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 refractionprofiles. Direct wave traveltimes are also included in each plot. Shots are located at eachend of the recording spreads and the maximum source-receiver offset distance is 50 m. Theprofile 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).(ms)(ms)(ms)(ms)N)-NC),000000000000000000000000Co qqCC)‘t.CDN) 0p,CDC; (tp,••—.‘bDc-f-p [0aCDLIIC’‘Cn‘—sCD:.:CO‘CDC‘—4CDCDCI) C H63The straight line arrival time curves plotted in each panel convey an impression ofa two dimensional earth model. The full three dlimensionality of the subsurface is onlyappreciated by comparing traveltime curves recorded along several separate azimuths. Notethat the reciprocal times (the shot-to-shot traveltimes) for forward and reverse arrivals in allpanels 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 traveltimecomputation method does not locate the critical offset distance, equation (3.32) is evaluatedover the full offset range covered by the receiver array. If traveltime analysis is concernedsolely with first arrivals, then these calculated nonphysical traveltimes do not pose anyproblems, because they are always associated with later arrivals. However, precritical offsetarrivals do have interpretive significance (Ackermann et al., 1986) and thus their inclusionin 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 timesobserved 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 dimensionalmodel consisting of a single layer overlying a halfspace, the intercept time is also independentof recording proffle azimuth (see equation (2.10)). Dependence on the azimuth angle ‘ onlyarises when multiple layers in three dimensions are analyzed. However, raypath reciprocityrequires that the propagation direction vectors satisfyPk(T! + ir)—qzk(1’), qk(P + ir) = Pik(’1’). (3.33)Substituting these results into the expression for the intercept time yieldsys, + 7r) bk(XS, ys, P). (3.34)64Hence, the two dimensional interpretive rule also holds for the particular three dimensionalmodel examined here (uniform velocity layers bounded by plane interfaces). However, itis not true that that intercept times recorded on all line profiles emanating from the sameshotpoint are identical. For realistic earth models (maximum dip 100) this variation inintercept 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 geometryThe final example examines head wave traveltimes recorded in a simulated VSP experiment. The computation procedure described in section 3.6 is readily generalized to asituation where the source and receiver are located on different interfaces of the model. Anazimuth angle function ‘1= g() can be calculated for a hypothetical source located oninterface j and a hypothetical receiver located on interface 1. In an offset VSP survey, thesource is located on the surface (j = 1) and a borehole geophone is lowered continuouslydown the well. Thus, an azimuth angle function is computed for all interfaces 1, includingthe critically refracting horizon (1 = k). These functions are then inverted at the knownsource-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. Thetraveltime to a geophone located within a layer is then calculated by linearly interpolatingtimes 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.65240— 220(I)200180240‘—s 220U)200180240r- 2201.i)2001800 160Fig. 3.7. Head wave arrival times recorded in a VSP configuration. Surface source is offsetfrom the well by 500 m to the north, east, and west in panels (a), (b), and (c), respectively.40 80 120depth (m)66The 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 dimensional formulae are needed to compute accurate arrival times. Figure 3.7 displays head wavetraveltime curves as a function of geophone depth. within a borehole positioned at the coordinate origin. Surface sources are offset from the well by 500 m to the north, east, andwest (panels a, b, and c, respectively). For this 2D model, a source offset 500 m to thesouth generates traveltimes that are identical to those plotted in panel a. The arrival timeof a given head wave decreases as the geophone is lowered in the well, because the receiverapproaches the critically refracting interface more closely. As the geophone passes throughan interface, the slope of the traveltime curve changes. Expressions for these slopes can beobtained from equations (3.17) or (3.18) above. Each curve terminates at the depth of thecritical 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 initialarrivals. Other waves that are neglected in the modeling procedure (e.g., direct waves) mayactually arrive first over certain ranges of receiver depth.3.8 ConclusionThe traveltime formulae derived in this chapter are useful for a variety of forward modeling applications involving head waves propagating in three dimensional, multilayered earthmodels. Obviously, construction of traveltime curves for a trial earth model can assist in theinterpretation of field recorded data. The equations also provide means for evaluating theimportance of three dimensional effects on head wave traveltimes in various seismologicalcontexts. For example, Merrick et al. (1978) investigate the hidden layer phenomenon usinga 2D earth model. Hunter and Puilan (1990), using a 1D model, compare the sensitivitiesof vertical and horizontal receiver arrays for discriminating layer velocities. Many investigators are concerned about the possibility that head waves might constitute first arrivals incrosswell transmission tomography experiments. All of these studies can benefit from a fullthree dimensional analysis.67Several straightforward extensions of the results described in this chapter can enhancethe utility of the formulae for forward modeling purposes. Inclusion of multiple reflectionswithin 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 appliedto calculate the particle displacement amplitude and waveform of a head wave in this threedimensional situation. Care must be exercised in treating shear wave propagation, becauseSV and SR modes are not globally decoupled in 3D. Richards et al. (1991) examine some ofthese 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 formula, inverse methods designed to recover the earth model parameters from measured dataare facilitated. Chapter 4 describes an inversion procedure that exploits the rapid forwardmodeling capability developed above. Other head wave traveltime inversion procedures areprobably possible. In particular, the ‘slope and intercept’ equation (3.32) may allow the extension of the two dimensional methods of Dooley (1952), Adachi (1954), and Johnson (1976)to three dimensional models. Although considerable effort has been devoted to discoveringthis generalization, success has not yet been achieved.68CHAPTER 4INVERSION OF HEAD WAVE TRAVELTIMES FORTHREE DIMENSIONAL PLANAR STRUCTURE4.1 IntroductionNumerous investigators have studied the inversion of reflection traveltimes for threedimensional 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 analogoussituation for refraction traveltime data has not been thoroughly examined. Kanasewich andChiu (198) present a method for jointly inverting reflection and refraction traveltimes torecover 3D structure. Forward calculation of traveltimes is achieved with the iterative raybending method of Chander (1977a).This chapter describes an algorithm for inverting head wave arrival times for three dimensional planar structure. The earth model is characterized by a stack of uniform velocitylayers bounded by plane, dipping interfaces. The inversion method is iterative; an initialestimate of the model parameters is refined until an acceptable match is obtained betweenobserved and predicted traveltimes. Rapid forward modeling of both traveltimes and trayeltime derivatives is achieved with the numerical technique developed in Chapter 3. A novelfeature of the inversion procedure is the inclusion of constraint information in the form of inequality relations satisfied by the model parameters. Often, a priori geological or geophysicalinformation is available to guide and constrain a traveltime inversion. This is particularlyuseful for the inversion of head wave arrival times because the problem can be very ifi-posedand admit numerous solutions.After a discussion of the mathematical basis of the inversion technique, the algorithmis tested on both synthetic and field acquired traveltime data. Inversion of refraction dataacquired in the Peace River Arch region of northern Alberta indicates that the algorithmcan be a valuable tool for the analysis of traveltimes recorded in a broadside configuration.694.2 Inversion mathematics4.2.1 General theoryLet the observed arrival times for the experiment be organized into an M-dimensionalcolumn vector t0. The earth model is characterized by a finite set of scalar parametersm, j = 1,2,. .. , N. These are organized into an N-dimensional column vector m. Predictedhead wave traveltimes generated by this model are designated by the M-dimensional vectortd(m). Then, a first order Taylor series expansion of the observed data about the particularmodel m yields the expressionA(m’) zm”’ = t(m), (4.1)where zt(m’) t0,3 — t,,.d(m) is the data discrepancy vector, m1 mt1 — m’ isthe parameter update vector, and the elements of the M x N sensitivity matrix A(m?z) aregiven by[A(m’)J = 2 (4.2)23 8mIn these and subsequent expressions, the superscript n denotes the iteration index.In many crustal seismic reflection and refraction experiments, the system (4.1) is overdetermined, underconstrained, and inconsistent. A popular solution technique for m’1 isthe damped least squares method (Braile, 1973; Kanasewich and Chiu, 1985; Chiu et al.,1986; Chin and Stewart, 1987; Phadke and Kanasewich, 1990). The updated parametervector is then obtained via m’1 = m’ + Zmfl+l. Iterations continue until an acceptablefit to the observed traveltime data is achieved. This strategy is termed ‘creeping’ (Scaleset at, 1990) because the final solution is obtained by the addition of (possibly many) smallperturbations to an initial guess.An alternative approach is used in this study to obtain the improved model parametervector. Substituting the definition of the parameter update vector into (4.1) and rearranging70terms givesAt(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) referto this technique as ‘jumping’ because, in the absence of additional constraints, the size ofthe model change between iterations m — JJ is not restricted to be small. Since theinversion is formulated in terms of the model itself, rather than a model perturbation, thejumping strategy facilitates the incorporation of constraint information into the algorithm.Constraints are mathematically expressed in the form of inequality relations satisfied by themodel parameter vector:m_ < m1 m, (4.4a)where the vectors m_ and m+ are lower and upper bounds on the model, respectively. Thesebounds arise from a priori geological or geophysical knowledge (or assumptions) about theearth 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 theinversion procedure. This is required for meaningful forward modeling of traveltimes. Theinequality bounds in (4.4a) can also be used to severely restrict (or even eliminate) thevariation 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 betweeniterations. Hence, if Sm is a vector of upper bounds on the parameter increments, theupdated parameter vector must satisfy the additional inequality constraintsm — Sm < m < m + Sm. (4.4b)These constraints fulfill the same regularizing role as the damping parameter in a dampedleast squares solution of the original equation (4.1). If a reasonably good initial estimate71m0 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 thanjumping 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. Arobust solution to this inconsistent system can be obtained by minimizing the l norm ofthe misfit. Linear programming provides a convenient solution method because the modelparameters are intrinsically non-negative and are constrained by inequalities (4.4a,b). Inorder to pose the problem in the context of linear programming, an M-dimensional residualvector 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. Theproblem now consists of determining the model parameter vector m1 and the residualvector r that simultaneously satisfy the inequality constraints (4.4a,b), the equality constraints (4.4c), and that minimize the 1 norm of the residualr = (4.4d)A standard linear programming routine is used here to solve the constrained optimizationproblem specified by equations (4.4a,b,c,d). After an improved model parameter vectorm’ is obtained, the lj norm of the misfit between observed and predicted traveltimestob3 — t,.d(m’’) is computed. Iterations cease when the misfit reaches some acceptablelevel, or exhibits negligable change on successive iterations. The elements of the residualvector 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 howclosely the linearized data equations have been fit.72Finally, the flexibility of the inversion algorithm is enhanced by including variable weighting 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, theseare restricted to be diagonal matrices. Thus, premultiplication of the sensitivity matrix byWd corresponds to row weighting and postmultiplication by W’ corresponds to columnweighting.The data weighting matrix can be used to emphasize those particular traveltimes judgedto be more significant for the inversion. The parameter weighting matrix serves to nondimensionaiize and normalize the elements of the model parameter vector m1. This is apractical concern in inversion algorithms where the model is characterized by parameterswith different physical dimensions and/or widely varying numerical magnitudes. The conditioning 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. Suitable units of measure are chosen for the various model parameters; the reciprocals of thesescalars form the diagonal elements of Wi,. The inequality constraints (4.4a,b) on the modelmust also be nondimensionalized in the same manner. If the weighted parameter vectorcalculated by the inversion algorithm is designated thL+l, then the physical parametersrequired for forward modeling of traveltimes are obtained via m’ = W;1 th’’1.4.2.2 Calculation of sensitivitiesA simple 3D earth model consisting of a single layer overlying a halfspace is characterizedby the five-element parameter vector m = [vi,v2,qf,O,hjTwhere v and v2 are P-wavevelocities of the two media, and 0 are interface orientation angles, and h is the verticaldepth to the refractor below the coordinate origin (Figure 2.1). A closed form expression forhead 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 straightforward73partial differentiation. However, a substantial simplification arises by redefining the setof 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 inthe slowness Si. In terms of these new parameters, head wave traveltime isT = 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). Differentiating 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 relationswith 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 rankdeficient. 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, for74example, if a one dimensional earth model is used for the initial parameter vector estimatem0.Calculation of head wave traveltime sensitivities for a multilayered earth model is morecomplicated because a closed form mathematical expression for the arrival time does notexist. However, Chapter 3 describes a rapid computing procedure for obtaining criticallyrefracted traveltimes for arbitrary recording geometries. A combination of analytical andnumerical techniques can then be used to calculate the sensitivities. In. this study, sourcesand 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 interfaceorientation angles, and K — 1 layer thicknesses). These parameters are organized into theN-dimensional column vectorm = [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 forlayers and interfaces is described at the beginning of Chapter 3.A ‘slope and intercept’ expression for the traveltime of a head wave critically refractedat interface k (2 k K) is given by equation (3.12). In the current notation, it isT = m(’P)X + b(x,y,’JI), (4.6a)where the slope m(P) and intercept time bfrs, ys, II’) arem(’TF) = s [cos iJ qlk,z + sin ‘I’ qik,y — tan q cos(P—0’) ik,z], (4.6b)75ys, = 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 propagationvectors Pk and characterizing the raypath depend implicitly on both the model m andthe recording azimuth ‘E’. Nevertheless, if m and are specified, then all of the propagationvectors 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 identicallyzero. Moreover, the sensitivity to layer thickness h1 for 1 < k is derived directly from theabove expression for intercept time: ÔT/0h1 = Si (pi — quz). The remaining sensitivitiesmust be evaluated by a finite-difference technique. Let dm1 be a model perturbation vectorwith zeros in all element positions except the th• Then, the partial derivative of head wavetraveltime with respect to parameter m is approximated by the forward finite-difference07’ T(m + dm,) — T(m) (47)—IIdmiIIUnit propagation vectors are generated for both the perturbed model m + dm, and theunperturbed model m via the procedure described in Chapter 3. These are then usedin formulae (4.6a,b,c) above to calculate traveltimes for each model. Finally, substitutingthese traveltimes into (4.7) yields the required derivative. The size of the model perturbationdm, is typically about 1% of the value of the associated model parameter. Although acentered finite-difference scheme would yield greater accuracy, the one-sided approximation76used here requires less computational effort. Only two traveltime calculations are necessaryfor each model parameter, instead of three.4.3 Synthetic examplesThe two examples discussed in this section are representative of a large number ofcomputational experiments conducted with the inversion algorithm. Both the single-layerand multiple-layer variants of the algorithm are examined. The first example demonstratesthat the single-layer version is capable of returning the correct solution under a variety ofoperating conditions. However, as indicated in the second example, the multilayer versionappears to require fairly restrictive constraints in order to iterate toward the correct model.4.3.1 Single layerThe earth model used for the first example consists of a single layer overlying a halfspaceand 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 equation (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 receiversare deployed around the perimeter of an equilateral triangle with sides 500 m long. Eachside contains 11 receivers separated by 50 m. All receivers record energy from a source thatis activated sequentially at the three vertices of the triangle. However, some source-receiveroffsets are less than the critical offset distance (equation (2.11)) and thus the receivers donot detect a head wave arrival. These fictitious times are excluded from the inversions. Atotal of M = 73 uniformly weighted traveltimes are used to recover the N = 5 earth modelparameters.77I100 + +0’ + ++++++ +++*+ + + + + + + + +*—200I I I I I—200 0 200200b++++100-+ 1+—200 0 200easting (m)Fig. 4.1. Plan views of two areal recording arrays. Sources are indicated by asterisks andreceivers by small crosses. (a) Triangular geometry. (b) Swath geometry. The strike andclip symbol in the upper right corner of each panel refers to a subsurface interface located100 m below the coordinate origin (large cross).78Numerical results from a few typical inversion runs are described here. The iterativeinversion procedure is initiated with a one dimensional model given byvi =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 terminatedwhen the relative change in the misfit is less than 1%. A wide variety of starting modelsyields essentially the same final solution, although the number of iterations required forconvergence varies. Also, the inversion is stable when the traveltimes are contaminated withsmall 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 startingmodel. 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 thenoise. Note that the overburden velocity vi has been correctly estimated from the refractiondata alone! This interesting (and unusual) result is consistent with the theoretical analysispresented in Chapter 2.These results suggest that the triangular recording array is a useful tool for determiningthree dimensional planar structure. This particular geometry combines an adequate distribution of offset and azimuth with three reciprocal time pairs. These are favorable attributesfor a successful inversion of refraction arrival times via the time term method (Scheideggerand 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 that79it yields an exact solution for the delay times at the three vertices. More recently, the triangular array has been used for deep crustal seismic exploration in Saskatchewan (Kanasewichand Chiu, 1985) and British Columbia (Zelt et aL, 1990).Other recording geometries are considerably less robust in detecting and resolving thethree dimensional dipping structure. Figure 4.lb depicts two parallel line arrays, separatedby 200 m, oriented along the strike direction of the subsurface refractor (Nw-SE). Eachspread contains 11 geophones (receiver interval 50 m) that record head wave arrivalsfrom a source located at the center of the opposite array. This broadside recording patternsimulates aspects of the ‘swath geometry’ commonly used for 3D seismic reflection surveys.Using the same starting model, inversion of the 22 broadside traveltimes yieldsvj = 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. Althoughthis model generates an exact fit to the data, the overburden velocity v and vertical depthh are incorrect. A similar effect is observed when the inversion is initiated from numerousdifferent starting models. This situation illustrates the classical tradeoff between overburdenvelocity and refractor depth in seismic refraction interpretation. If additional a priori dataare introduced into the inverse problem, then a correct solution is possible. For example, theinterface depth may be known from a borehole drilled at the coordinate origin. Constrainingthe 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 theinequalities 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.80Evidently, the broadside geometry illustrated in Figure 4.lb allows many solutions to thisnonlinear traveltime inverse problem. However, if the inline head wave arrivals recorded byeach spread are also included in the inversion (M = 32) then the correct model is recoveredin 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 refractingmedium. Equation (2.13) and Figure 2.4 indicate that the measured apparent velocity equalsthe true velocity in this situation. As expected, the inversion is degraded when the exacttraveltimes are contaminated with uniformly distributed random errors on ±4 ms. Themodel= 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 layersThe earth model used for the second example consists of two layers overlying a halfspaceand is defined by the parameter valuesv = 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 theinversion. Hence, there are only N = 9 parameters to estimate. Head wave traveltimes for81each critically refracting horizon are generated by evaluating equation (4.6a). The recordinggeometry used is the triangular array displayed in Figure 4.la. Also, precritical offset arrivalsare included in the traveltime dataset (in an actual field experiment, these times could beestimated by extrapolation or phantoming). Hence, there are a total of M = 174 arrivaltimes input to the inversion procedure.If restrictive constraints are imposed on the model parameters, then it is possible to recover the correct solution with the multilayer inversion algorithm. For example, constrainingthe 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 areplaced 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 relativelyunconstrained 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 linearized data equations, has many local minima that preclude convergence to the desiredglobal minimum.4.4 Field data example441 Peace River Arch broadside dataDeep seismic refraction data were acquired in the Peace River Arch (PRA) region ofnorthern Alberta in 1985. In addition to four inline profiles, two broadside profiles wererecorded. Figure 4.2 illustrates the source-receiver geometry for these two proffles. ShotA4 in the west is recorded by the north-south trending line A, and shot B4 in the north isrecorded 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 wavesthat are critically or near-critically refracted at the Moho (Zelt, 1989). Hence, inversionof the first break traveltimes can provide an estimate of the regional depth and dip of theMoho beneath the Peace River Arch.82-cCC*Shot 84300-200-100 -0—100—200Line BShot A4 Line A—200 —100 0 100 200easting (km)Fig. 4.2. Broadside recording geometry for the Peace River Arch seismic experiment.Sources are indicated by asterisks and receivers by small crosses.83In 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 pictureof the crust and Moho in the Peace River Arch region. Rather, the intent is to recover a largescale 3D model that can be subsequently refined by other traveltime interpretation/inversionmethods. Initial estimates of the model parameters are obtained from crustal sections alonglines A and B given by Zelt (1989, p. 103 and 112). These sections were derived byinterpreting the inline refraction data of the PRA experiment via a trial-and-error forwardmodeling 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 estimatedpicking error is relatively large (±50 ms).442 Static correctionsThe effect of variable near surface structure on the head wave arrival times can be reducedbe applying static corrections to the first break picks. Static corrections are commonlyapplied to seismic reflection data for this same purpose. However, there is an importantdistinction between the corrections for the two types of data. In the refraction case, thestatic 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. Inthe reflection case, the correction pertains to the vertical traveltime through the actual andreplacement media. Statics application is an important preprocessing step for the PRAbroadside data because the subsequent inversion assumes a very simple earth model. Ineffect, an attempt is made to ‘make the data fit the model’ more closely. The followingdevelopment assumes that near surface velocity information is available. For the PRAexperiment, 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 approximatingthe local near surface velocity structure by a one dimensional stack of layers. Theth layeris bounded from above and below by plane, horizontal interfaces at depths z, and zj,84respectively. The P-wave velocity of the layer depends linearly on depth z and is given byv(z) v(O) + a,z. In this situation, the raypath associated with a particular ray parameterp is a circular arc. The upward propagation time across the layer isit(p) = -- cosh (_2_) — cosh1 (—_) , (4.8)a1 pv1 pvwhere v v(z) = v.(O) + az and v vj(zj+i) = v(O) + a.jzj+1. Symbols v and v denotethe velocities at the top and bottom of the i layer, respectively. The horizontal distanceaccumulated by this wave traversing layer i isIXj(p) =—-- [‘i - (pvj)2 - - (PvD2]. (4.9)If the wave is critically refracted on a horizontal interface at greater depth, then the rayparameter p equals 1/va, where v is the critical refraction velocity. The angle i(z) thatthe 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 awhere ?7 = sin(v/v) and = sin1(v/v). 77 and 77 are the incident angles of thecircular raypath segment at the top and bottom of layer i, respectively.The refraction delay time contribution of the th layer is defined asJt —85Substituting from equations (4.10) yields= ---ln , (4.11)a1where 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 velocitygradient of the layer vanishes (a —* 0). In this situation, the raypath becomes a straightline 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 thicknessh and uniform velocity v(0).The one-way refraction delay time associated with a stack of m horizontal layers isobtained 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 PRA86example), then Q = Q+i and expression (4.12) simplifies to= E--1n . (4.13)a Q+iThe 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 delaytime associated with a uniform replacement medium of equal thickness. If the replacementvelocity 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 ofupward wave propagation through a set of horizontal layers, it also applies to the downwardpropagating portion of the total raypath. Hence, it may be used for the computation ofeither a source site static or a receiver site static.Near-surface velocity information along lines A and B of the PRA experiment are obtained from models developed by Zelt (1989, p. 40) from well logs. There are eight wellsalong line A and seven wells along line B. At each weilsite, the velocity model is approximated by a one dimensional stack of layers with linear velocities, and a static correctionis calculated via formula (4.14). The critical and replacement velocities are assumed to be8.25 km/s and v,. = 6.6 km/s, respectively. These weilsite static corrections are theninterpolated/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 weilsitesalong a proffle using the York algorithm (York, 1966, 1967, 1969; Williamson, 1968). Thisis the proper line fitting technique to use in this situation because both the northing andeasting coordinates are subject to positioning errors. The York algorithm also automaticallyreturns the coordinates of the perpendicular projections of the points to the fitted straight87line. Finally, the weilsite statics are linearly interpolated/extrapolated in these projected coordinates to the receiver sites. Figure 4.3a displays the computed receiver static correctionsfor both lines; squares denote the weilsite statics calculated from expression (4.14). Themagnitudes of the statics are approximately the same on both lines ( 300 ms), althoughthere 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 readilyestimated for this shot from the nearest computed receiver correction (—290 ms). A sourcestatic for the isolated shot B4 must be assumed; the value adopted here is —300 ms.4.4.3. Inversion resultsInitial estimates of the model parameters are inferred from the interpreted crustal sections in Zelt (1989, p. 103 and 112). These starting parameter values arev = 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 modelparameters: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 parametersSj (slowness) and i (critical angle) that are used by the inversion algorithm. After fiveiterations, the following model is returned:= 6.21 km/s, v2 = 8.50 km/s, = 2.75°, 8 = 299.74°, h = 42 km.EE88o ‘100 200 300—0.20 -0 100 200 300—0.20 L:0 100 200 3000 100 200 300inflne distance (km)Fig. 4.3. Receiver static functions for lines A and B of the PRA experiment. Squares referto well locations where near surface velocity logs are available. (a) Statics calculated withv = 8.25 km/s and yr = 6.6 km/s. (b) Statics calculated with v = 8.5 km/s and v1. = 6.2km/s.89The initial traveltime misfit of 658 ms is reduced to 172 ms; iterations cease when the relativechange in the misfit is less than 1%. The velocities vi and V2 of this solution differ fromthe 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 appliedagain to the picked arrival times. Figure 4.3b ifiustrates this second set of receiver staticcorrections. New source statics are assumed to be —210 ms for shot A4 and —250 ms forshot 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 differentfrom 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 character 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 subsurface structural complications. However, a large component of the total traveltime misfitmust be attributed to structure or velocity variations that are not modeled in the inversionprocedure. This misfit is too large to be accounted for by random picking errors alone. Forexample, the predicted times are systematically greater than the observed times throughoutthe central portion of line B. This suggests that the Moho north of line B is not adequatelyrepresented by a plane interface bounded by uniform velocities.Vertical depths to the Moho calculated from the inversion results are posted on a planview 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 dipof 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 recoveredby 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.90Line A0048464442403850484644424038I I I I I0 100 200 300Line BI I I I I I I0 100 200 30091Fig. 4.5. Vertical depths (km) to the Moho in the PRA region inferred fom the inversionresults.300200-100 -0-—100 -—2000’-c-I-?0C42.*:34.642.047.237.747,8I I I—200 —100 0 100 200easting (km)92precise 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, Zeltinferred a P velocity of 8.4 km/s along the northern half of line A and the eastern quarterof line B.An independent check on the validity of the 3D inversion is provided by a methodproposed by Zelt (1989). Using a simple two dimensional earth model, each arrival timerecorded along a broadside proffle can be inverted for an estimate of Moho depth beneaththe associated receiver site. A depth proffle for the Moho can then be constructed byplotting these depth estimates side-by-side. Zelt (1989) accomplished this inversion usingray tracing techniques. However, as the following derivation demonstrates, a straightforwardmathematical solution to the problem is possible. In the two dimensional situation, headwave traveltime isT(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 thelayer 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 yieldsT(X) = sin(i+’p) + 2h(O) C:: ‘p cos c + 2 sin’p cosh(O) is the layer thickness at the coordinate origin. This relation is rewritten as a quadraticform in cosAcos2’p + Bcos’p + C = 0. (4.16)93The three coefficients depend on earth model parameters {vl,v2, h(O)} and measureddata {xs,X,T}:A = 1 +4h(o)cosi [ni + [1 + (i + )jJv1T F 2h(O) 1B = —2—i— Lsmn z + x cos2viTi 2x 2c=—1+_ cosic.If values for v, V2, and h(O) are known (or assumed), then these coefficients can be evaluatednumerically. Solution of equation (4.16) for cos is via the quadratic formula:-B+/B24ACcos(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 ambiguityis easily resolved by ensuring that the traveltime predicted by formula (4.15) agrees withthe known traveltime T. Finally, once the dip angle is determined, the vertical depth to therefractor 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 inferredfrom the parameters recovered by the 3D inversion. Vertical depths beneath the recordingstations on each line are plotted. The coefficients in the quadratic (4.16) are evaluated withthe 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 variationsin 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 completely94different techniques, especially along line A. The larger departure of the two depth curvesalong 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 reasonable depth estimates, it requires the assumption of numerical values for the three unknownmodel parameters {v, V2, h(O)}. In contrast, the 3D inversion technique yields simultaneousestimates of all of the relevant earth model parameters. Moreover, since it accomplishes ajoint inversion of all of the error contaminated traveltime data, it is more robust than the2D method.4L1neAO 100 200 3000 100 200 300in’ine distance (km)Fig. 4.6. Comparison of vertical depths to the Moho beneath the receivers of lines A and Bcalculated by two different methods. Smooth curves in each panel are depths inferred from3D inversion results. Jagged curves are obtained from a simple 2D inversion method.4.5 ConclusionThe inversion of head wave arrival times for three dimensional planar structure is posedas a constrained, parameter optimization problem. The iterative inversion algorithm de95scribed here has the abffity to converge to a realistic solution provided that (i) the data acquisition 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 geometries for detecting and resolving three dimensional dipping structure. Also, investigation ofboth synthetic and field recorded datasets indicates that it can be successfully applied to theinversion of broadside refraction data. Currently, there is a lack of techniques for effectiveinterpretation of such data. Although the multiple-layer version of the algorithm exhibits agreater tendency to converge to an erroneous result, several successful inversions have beenachieved with the inclusion of sufficient constraints. Generalization of these algorithms tothree dimensional recording geometries would be a useful extension. As indicated in Chapter3, the forward modeling procedure can easily accommodate buried sources and/or receivers.Moreover, well log information regarding interface depths and layer velocities would providethe necessary constraint information for the multilayer algorithm.96CHAPTER 5POINT DEPTH AND DIP ESTIMATES FROMREFRACTION TRAVELTIMES5.1 IntroductionThe Generalized Reciprocal Method (GRM) is a technique for delineating an undulatingsubsurface interface from refraction traveltime data recorded by inline forward and reverseproffles. After the refracted arrivals from a particular horizon are identified, simple linearcombinations of the picked times are used to form two analysis functions: (i) the velocityanalysis function, and (ii) the generalized time-depth function. Each is plotted with respectto receiver position along the profile line. The slope of the velocity analysis function is ameasure of the P-wave velocity of the refracting layer. The generalized time-depth functionis 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. Conversionof the time-depth function to actual depth values requires velocity information for both theoverburden and the refracting layer.The GRM was developed by Palmer (1980, 1981) as an extension of the conventionalreciprocal time method for interpreting refraction arrival times (Hawkins, 1961). Althoughit has been successfully applied to the problem of mapping undulating refractors, the mathematical formulation of the method is based on a two dimensional earth model with planeinterfaces. 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 thicknesses is “not suitable” for the derivation of the ORM analysis functions. This statement isincorrect. All of the GRM analysis tools can be rigorously derived from the two dimensionalversion of equation (3.23) for head wave traveltime. This development is the subject of thefirst portion of this chapter. The use of vertical layer thicknesses results in substantiallysimplified derivations. Hence, this reformulation of the theoretical basis of the GRM is quite97useful for educational purposes. However, a practical benefit related to refractor depth profile construction is also realized. The equations developed here allow point depth estimatesof the interface to be calculated from the generalized time-depth function. Computationof a refractor depth profile then reduces to an interpolation problem (treated in Chapter6). In contrast, Palmer’s time-depth function yields a circular locus of possible refractorpositions. Depth profile calculation then involves the more complicated task of constructingan envelope to a set of circular arcs.As indicated above, the GRM is founded upon a two dimensional earth model with planeinterfaces. A traveltime inversion technique that explicitly incorporates nonpiane refractinghorizons into the model should yield more accurate results. Hence, the second portion of.this chapter proposes a 2D head wave traveltime interpretation method tentatively namedcritical offset refraction profiling. This inversion technique accommodates undulations inboth the refracting interface and the surface, as well as horizontal variations in the velocityof the refracting layer. Detection and definition of such structural and velocity variationsare subjects of current research (Palmer, 1991). Critical offset refraction profiling yields amathematically exact solution to the 2D seismic refraction problem; point values of interfacedepth, interface dip, and refractor velocity are obtained. A continuous refractor depth profilecan then constructed by interpolation. The dip estimates provide additional constraints onthe interpolant.5.2 Specialization to a 2D modelPrior 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, theorientation of interface i is described by a dip angle qfj (0 ç& <ir/2) and an azimuth angleS (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 isn1 = (sin gj cos O)i + (sin çt sin 02)j + (cos (5.la.)98and the vertical depth of the interface isz(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) ish(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 restrictedto the two values 0 or ir. In this case, the above equations reduce ton = (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) areindependent of the y-coordinate. For subsequent analysis of the two dimensional refractionproblem, it is convenient to repararneterize the earth model in terms of interface dip anglesthat may be positive or negative. Hence, let the symbol Oj refer to an interface dip anglein a two dimensional layered earth. cpj is an acute angle (0 I < 7r/2) measured withrespect 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. Rewritingthis relation as = —çoj (cos 9i) and then substituting into equations (5.2a,b,c) gives thefamiliar expressions99= (— sin p)i + (cos y)k, (5.3a)z1(z) = zj(O) + x tanyj, (5.3b)h(x) = h(O) + sin(cpjl— (53)COS COS çoThe head wave traveltime formula relevant to GRM analysis is equation (3.23). Thisis specialized to a two dimensional model by setting ys = YR = 0, n = — 5111’pk, andk,z = cos c°k. Equation (3.23) becomesT,(xs, XR) =—+—COS Yk + Pkk,z SilL Yk) (5.4)vkcoscokA further simplification is obtained by expressing the components of the unit propagationvectors Pik and q in terms of orientation angles. In the three dimensional situation, eachof 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/3jare poiar angles measured from the +z axis (0 ak, I3ik ir), and and8ik are azimuthal angles measured from the +x axis (0 7ik, < 2ir). If the recordingprolile is oriented perpendicular to the strike of the subsurface horizons, then afl propagationvectors are confined to the xz plane. The azimuth angles equal 0 or r, and expressions(5.5a,b) reduce toPik = (+SiflUk)i+(C0Sak)k, q:k (+sin/3k)i+(cos/3k) . (5.6a,b)100The propagation direction of the critically refracted raypath segment is obtained by straightforward 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 resulth(xS)cosak — h(xR)cos/3k IZR — sITk(XS,XR) = + . (5.7)i=1Vk COS PkEquation (5.7) is a novel expression for surface-to-surface head wave traveltime thatforms the basis for an alternative development of the GRM. It differs from the analogousequation given by Palmer (1980, p. 5) in several important ways: (i) layers are characterizedby vertical thicknesses below source and receiver, (ii) raypath orientation angles are measuredwith respect to the vertical, (iii) the coefficient multiplying the source-receiver offset distancedepends only on critical refractor quantities (i.e., velocity vk and dip and (iv) the earth’ssurface may be nonhorizontai. These attributes facilitate a straightforward derivation of thevarious GRM analysis tools.The raypath angles in equation (5.7) depend on the recording proffle azimuth , whichin turn is restricted to the two values 0 or r. However, raypath reciprocity requires thatak(7r) = —I3k(0) and f?(ir) — ak(0). These expressions imply thatcosak(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) satisfiessource-receiver reciprocity: Tk(XS,XR) =1015.3 Derivation of GRM parametersThe two dimensional earth model used to derive the GRM analysis functions consistsof uniform velocity layers bounded by plane, dipping interfaces. Figure 5.1 depicts therelevant critically refracted raypaths. Although the subsequent mathematical developmentassumes 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 andZs3, respectively; 1? denotes a typical receiver with coordinate ZR. The velocity analysisfunction and the generaiized time-depth function are constructed from the traveltimes offorward and reverse propagating head waves that arrive at the positions ZR + L/2 andZR — L/2, respectively. L is an adjustable separation distance. If £ = 0, then arrivals atthe single receiver at ZR are examined, and the GRM reduces to the conventional reciprocalmethod. GRM analysis seeks to determine the particular separation £ such that the upgoingsegments 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-depthfunction 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 describedby Palmer (1980, 1981, 1990, 1991) and are not investigated here.5.3.1 Velocity analysis functionThe velocity analysis function for layer k is defined asT(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).102Fig. 5.1. Critically refracted raypaths used to derive the GRM analysis functions. R denotesa receiver located between the two surface sources S1 and S2.Substituting the head wave traveltime formula (5.7) into this definition yieldsk—ivii=1—h(:s2)[cosaik(ir)+cos/3ik(O)](5.10)0 XS1LXRV2+h(xj — L/2) COs/3jk(lr) — h(zR + L/2) COS/3k(0)+—xs12v t7k cos çoi=1where the dependence of the raypath angles on recording profile azimuth is explicit. Usingequation (5.3c), the layer thicknesses at the horizontal coordinates XR ± L/2 are written ash(xR ± L/2) = h(xs1)+ (XR — Xs1 ± L/2)sin(yj1—cos Pi+1 cos c°103Substituting this expression and relations (5.8a,b) into (5.10) and reducing yield the resultT(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, thevelocity analysis function is a linear function of the distance (ZR — xs1). The slope of thisline is used to estimate the velocity of the critically refracting layer. The intercept has asimple interpretation that is given in section 5.3.5.5.3.2 Apparent horizontal velocityThe apparent horizontal velocity of layer k is defined as1 = 8T(z; L) (5 12)Va ÔZRUsing 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 2Thus, the apparent horizontal velocity of a refracting layer is influenced by all velocitiesand dips in the overburden, as well as by the dip of the critically refracting horizon. However, for dip angles commonly encountered in refraction exploration, formula (5.12) yields areasonably 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) shallowsubsurface model is quite useful for this purpose. Since the interfaces have unusually large104dips, this model provides a strong test of the utility of (5.12) for velocity estimation. Furthermore, 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 parameterszi(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 thespecified geometric parameters). In order to evaluate equation (5.13), the raypath anglesak(O) and/3k(0) must be determined. Application of Snell’s law of refraction at interface iyields the two recursion relations‘i—1,k = sin1 [!z1 sin(ak + co)J — 13i—1,lc = ir — sin[Vi_1sin(/3k+ oi)j—Yi.(5.14)Starting with akk(0) + Yk =/3kk(O) + y = ir/2, these relations are evaluated successivelyfori=k,k—1,...,3,2toobtaina2 = 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) yieldsk=2: Va1827.77m/s, k=3: Va=2578.59m/s, k=4: Va4209.65m/s.105These 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 approximation 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, animproved estimate of refractor velocity can be obtained if the dip ço of the critically refracting interface is known. This estimate is derived as follows. Assume that the dip angles of allinterfaces overlying the critically refracting horizon vanish (pj = 0 for i = 1,2,. . . , k — 1).Then, equation (5.13) becomes1 1 — tan yj [cos ak—1,k(O) + COS /3k—i ,k(0)Va vkcoscpk Vk_i L 2Geometric analysis reveals that the raypath angles within the layer immediately aboveinterface 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 relationsinto the above expression for Va and reducing yieldVavk (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 thecritical horizon are small, rather than identically zero. A dip corrected apparent horizontalvelocity can be defined as V cos Pk Using the above model parameters, the followingvalues 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 are0%, +2.8%, and +1.7% for layers 2, 3, and 4, respectively. Note that VL calculated forlayer 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.106As indicated previously, the condition for exact validity of equation (5.15) is the dipangles 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 interfaces (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 velocities v = 1000 + 100 i rn/s and dips ço i°, ( i = 1,2,.. . ,50). All differential dip anglesequal 10 and thus are small. A comparison between the apparent velocities and the truelayer velocity is illustrated in Figure 5.2. In particular, both Va and V depart appreciablyirom the correct velocity value as the layer index increases, i.e., as the dips of the boundinginterfaces become large. This behaviour confirms that the approximation (5.15) requiressmall absolute dip angles, rather than small differential dip angles.Fig. 5.2. Comparison of apparent velocity Va and dip-corrected apparent velocity V withthe 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 setof parallel, but dipping layers. Assume that all interfaces (including surface and criticalrefractor) have the same dip angle (ço = ço, for i = 1,2,.. . , k). Obviously, the differentialU)10 20 30 40 50layer number107dip 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 fromthe 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 specializingequation (3.12) to a two dimensional earth model. Using this equation in the definition ofthe velocity analysis function yields1—8T — mk(O)+mk(lr) 517VaOR 2 •Thus, the apparent horizontal slowness is merely the arithmetic mean of the slopes ofthe forward and reverse traveltime curves. This implies that Va is the harmonic mean ofthe velocities measured on the forward and reverse spreads. Defining vf 1/m(O) andVr 1/mk(Tr), then1 11 1 VfV,.—= ——+ — , or V0 = 2V0 2 Vf Vr Vf+VrLankston (1990) and Palmer (1990) state that the harmonic mean of v1 and v7 equals thetrue layer velocity. The present analysis reveals that claim to be erroneous. Rather, theharmonic mean V0 is an approximation to the true velocity vk.1085.3.3 Generalized time-depthThe generalized time-depth function for interface k is defined asTd(XR; L) [Tk(zs ,ZR + £/2) + Tk(ZS2,ZR — L/2) — , xs2)— ]. (5.18)Substituting from equation (5.7) for the head wave traveltimes givesTd(ZR;L)[cosak(0)_cosik(0)]I, —1L j 1 ç- 1 s1n(+1 — yi) cosak(0) + cos8k(O) 1+2jvlccoScok 4vi CoSyjcoScpj 2 VaHowever, from equation (5.13) for the apparent horizontal velocity, the term in braces isidentically zero. The above expression reduces toTd(ZR;L)=h(x1) [cosaik;cosf3ilC] (5.19)where the explicit dependence of the ray angles on recording proffle azimuth has been suppressed because, from relations (5.8a,b), the difference cos k(I1)_cos/3(1P) is independentof azimuth. The time-depth of a critically refracting interface at coordinate ZR is related tothe vertical thicknesses of all overlying layers at the same position. Plotted time-depth sections give a qualitative indication of refractor structure, in a manner analogous to reflectiontime sections. Conversion of the time-depth to actual depth requires velocity information,and is briefly discussed in the next section.109The slope/intercept formula for head wave traveltime also yields an alternative expression for the time-depth function. Substituting equation (5.16) into the definition of thegeneralized time-depth and reducing givesTd(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 formulamk(0)(x52 — xs1) + bk(x51) = mk(lr)(x52— z51) + bk(x52).Hencebk(zS1) — bk(xS2)mk(7r)—mk(O) =xS2 —Finally, substituting this result into the above expression for the time-depth function yieldsTd(cR;L) = bk(z51) Zs2 XR . bk(a52) XRS1 (5.20)2 x2—2 Xs2 — Xs1The 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 classicalrefraction traveltime analysis (Wyrobek, 1956).1105.3. Depth conversion factorExpression (5.19) for the generalized time-depth may be rewritten asTd(XR;L) =h(x) (5.21)Vjkwhere the depth conversion factor is defined asVji,2v (5.22a)COS ajk — COSI3ikThe depth conversion factoi has physical dimension of velocity. Since the raypath angles exileand /3ik are usually not determined in a GRM interpretation, Vk cannot be evaluated numerically. However, various approximations to the actual depth conversion factor are useful. Ifthe cosines of the raypath angles are replaced with the expressions appropriate for horizontallayering (cos exile — cos = — (Vj/Vk)2), then equation (5.22a) is approximated byVjVk=_______.(5.22b)Alternately, the true layer velocities in (5.22b) can be replaced by apparent horizontal velocities to obtain-ai ale=__, (5.22c)i/VVawhere Vaj denotes the apparent horizontal velocity of layer i. Evaluating these expressionswith the parameters defining the shallow four-layer model yields the following numericalvalues (in m/s):111i k V Vik1 2 1221.23 1202.68 1194.661 3 1087.69 1091.09 1084.912 3 2603.66 2593.76 2591.171 4 1041.93 1032.80 1029.472 4 2013.08 2015.61 2029.003 4 3315.54 3202.56 3262.22The 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 observedtime-depths via expression (5.21). For k = 2, this equation yieldshl(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 thicknesses: zkfrR) = zl(XR) + hj(XR) for k 2. Note that the depth zk(XR) is anunambiguously defined function of the receiver coordinate ZR. In contrast, the ‘depth’ calculated via Palmer’s (1980) time-depth function is a radius of a circular arc centered at thereceiver position.Continuing with the illustrative example, numerical values of the time-depth for thethree subsurface interfaces can be obtained by evaluating expression (5.20). Sources arelocated at the horizontal coordinates zs1 = —100 m and xs2 = +60 m, and the receiver ispositioned at the coordinate origin. The results arek = 2: Td(0;L) = 16.63ms, k = 3: Td(0;L) = 33.74ms, k = 4: Td(0;L) = 45.23ms.112Finally, vaiues of interface depths at ZR = 0 are calculated using the exact and approximatedepth 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 obtained by using the approximate depth conversion factors.5.3.5 Time-depth near a shotpointThe intercept on the time axis of the velocity analysis function (5.11) has a simple, butnevertheless useful, interpretation. From formula (5.3c), the vertical thickness of the thlayer at the coordinate xs1 + £/2 can be written ashj(z1 + L/2) = h(xs1)+ sin(co:+1 —2 cos Pi+1 cos çoSubstituting this result into expression (5.11) for the velocity analysis function yieldsT1,(XR;L) ‘ h(xs1+L/2) [cosoik — COS/3jk]+ZR —ZS1However, from equation (5.19), the summation term is recognized as the generalized timedepth at the coordinate xs1 + L/2. ThusT(xR; L) = Td(ZS1 + L/2; L) + ZR ZS1 (5.24)113The 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 atool for estimating refractor depths near a source, where head waves are either nonexistentor difficult to pick accurately.5.4 Critical offset refraction profiling5.4.1 Earth modelThe two dimensional earth model considered in this section is ifiustrated in Figure 5.3. Asingle layer with uniform velocity v is underlain by a medium where velocity v2(z) dependssolely on the horizontal coordinate. Surface and subsurface horizons are undulating and aredefined by the depth functions zl(x) and z2(x), respectively. The dip angles of these twointerfaces are given by—1 Fdzi(z) dz2(x)= tan dx S02(Z)= tan dxThese relations define the sign convention associated with interface dip.Forward and reverse traveltime data axe recorded by receivers distributed along thesurface. The observed traveltime curves are designated Tf(z) and Tr(x), respectively. Theslopes of the two traveltime curves aredTf(z) — dTr()mf(x) dx m(x) = — dxNote that a negative sign is included in the definition ofm7(x). From the measured slopes ofthe 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 andreverse profile raypaths, respectively. Then, straightforward geometric analysis yields theexpressions114Fig. 5.3. A 2D earth model with undulating surface and refractor topography. Forward andreverse raypa.ths exiting from the same point C on the critical horizon intersect the surfaceat 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 previouslyestablished sign conventions for surface dip and traveltime slope determine the positivesense of these two incidence angles. In particular, jLf(c) is positive if the angle opens inthe clockwise sense. The converse holds for ILr(x); it is positive if the angle opens in thecounterclockwise sense (see Figure 5.4).0Lc(xM)XMziABviV2(X)115The measured traveltime data pertain to waves that are refracted at the subsurfacehorizon z2(r). In general, the raypaths of these waves may undergo noncritical refractionsat the interface. However, the inversion procedure discussed here adopts the following fundamental assumption: the critically refracted portion of the total raypath is coincident withthe undulating subsurface interface. Although this is an approximation, it is sufficientlyaccurate if the structure on the interface is not too severe. In the presence of syncinalstructure, the approximation is reasonable because there is a tendency for the diffractedraypath 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. Infact, it underlies the popular interpretive rule known as the ‘law of parallelism’ of refractiontraveltime 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 thesubsurface interface. If i2(z) sin1 [vl/v2(x)] is a variable critical angle defined alonginterface 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 designatedA 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 criticalrefraction can be observed. Obviously, the critical offset varies along the line of proffie dueto 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 offsetfunction L(xM) can be defined at every midpoint position on the surface. The inversiontechnique outlined below is based on the fact that an exact solution to the refraction problemis possible if this function is known. A procedure for estimating the critical offset functionfrom the measured traveltime data is discussed in a subsequent section.5.4.2 Inversion methodIf 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. Figure1165.4 indicates that— XB—XCtan ttf(XB) —zC — Z5XC — X4tan1t(xA) —zC — ZSolving these equations for the two unknown coordinates (XC, zc) yieldsXC =zC =(XB — ZA) COS /Lf(XB) cos liTsin (ti.e() +(5.25a)(5.251,)These relations establish the position of the point of critical refraction C in terms of knownquantities. Furthermore, since the sum of the angles of triangle ABC equals ir radians, thecritical angle at C is given bypf(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 ortriangle BCE in Figure 5.4. The result isS02(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)117Fig. 5.4. An expanded view of triangle ABC in Fig. 5.3. Angles j and p describe theorientation of raypath segments CB and UA relative to vertical, respectively. The criticaland 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 specifiedseismic refraction problem. The formulae are not restricted to the particular geometricsituation depicted in Figure 5.4, but remain valid if the refracting point C does not residebetween 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 refraction solution for an earth model with plane interfaces and uniform velocities. In thatcase, the forward and reverse traveltime curves are straight lines with slopes that are independent of the horizontal coordinate x. Hence, the forward and reverse angles of incidenceare also constant: ILf(x)=and p7(x) = pr. Equations (5.26) and (5.27) reduce to thefamiliar expressions i2=(ELf + p1)/2 and Y2 = (jLf — ,u)/2. There are no analogues toA118equations (5.25) in the classical solution, because vertical depths to the refracting interfaceare determined only at the source locations.A set of point estimates of the position of the critically refracting horizon is obtainedby 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. Additionalconstraints on. the interpolant are provided by the interface dip estimates at each refracting point C. Assuming an error-free inversion, the constructed depth proffle should passthrough the coordinates (ac, zc) and have slope tan cp2(xc) at this point. A procedure forconstructing a smooth interpolation of a set of point estimates of depth and dip is presentedin Chapter 6.5.4. 3. Critical offset determinationThe inversion method outlined above requires the critical offset distance L(xM) tobe 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 ona line. This section discusses an extension of their method to arbitrary locations along areversed refraction proffle. Moreover, an analytical, rather than computational solution tothe problem is given. Consider two surface locations with a common midpoint coordinateXM and separated by the horizontal distance L. A candidate critical refraction raypathcan be constructed from these two points by projecting the incoming raypath segmentsdownward until they intersect (Figure 5.5). The projection directions are determined fromthe known incidence angles f(rM + L/2) and (ZM — L/2). The propagation time alongthis hypothetical raypath is called a ‘predicted’ traveltime and can be derived via geometricanalysis. It is119T7,,.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 determinedfrom the measured forward and reverse traveltime curves as follows:XMA \1/BZ2(X)Fig. 5.5. Determination of the critical offset distance L(xM) at the midpoint coordinateZM. Raypath segments projected downward from A and B on the surface intersect at pointC on the refracting interface.T0b3(XM; L) = Tf(XM + L/2) + Tr(ZM — L/2) — TR, (5.29)120where TR is the reciprocal time: TR = Tf(xS2)= Tr(xj). For a fixed midpoint coordinateXM, expressions (5.28) and (5.29) are evaluated over a range of values of the separationdistance 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 commontime calculated at this separation distance is the critical offset traveltime Tc(XM) associatedwith the midpoint position. Application of this analysis technique at various midpointlocations along the profile line yields several point estimates of the critical offset functionL(x).5.5 ConclusionThe theoretical foundation of the generalized reciprocal method of refraction traveltimeinversion has been reformulated in terms of layers characterized by vertical, rather thannormal, thicknesses. This reformulation results in a simplification of the derivations of thetwo GRM analysis functions, and thus has educational value. However, it also allows pointestimates of interface depth to be calculated from the measured traveltimes. A refractordepth profile can then be calculated by interpolating the depth points. Although the GRMis based on a simple 2D model with plane interfaces, interpreters apply the technique inmore realistic environments with nonpiane refracting horizons. The results that are obtainedshould then be viewed as approximate.In contrast, the proposed method of critical offset refraction profiling explicitly incorporates undulating interfaces and variable refractor velocities into the model. A mathematically exact solution to the traveltime inverse problem indicates that point estimates ofinterface depth, interface dip, and refractor velocity can all be obtained from the observedarrival times. These point estimates are positioned properly in space; there is no need for asubsequent migration of the solution. In common with many other refraction interpretationmethods, independent knowledge of the overburden velocity v is required. Although thederivations 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 triangleABC in Figure 5.3. Lateral variations with a characteristic scale greater than the horizontali2 1dimension of the triangle are allowed. The primary approximation adopted by this techniquerelates to the refracted raypath: it remains coincident with the undulating interface. Theimplications of this approximation for the inversion should be investigated more thoroughly.Finally, extension of the method to multilayered earth models represents another avenueof developement. Using the present algorithm, this could be achieved via a layer strippingapproach.122CHAPTER 6CONSTRUCTION OF A SMOOTH REFRACTOR DEPTH PROFILE6.1 IntroductionLinear inverse theory provides techniques for constructing acceptable models that areconsistent with the available data. However, if acceptabffity is judged by the data misfitalone, then there are usually infinitely many equally valid solutions. A specific solution canbe obtained only by extremizing some functional of the model, subject to a requirementthat the final model is in satisfactory agreement with the observations. Since differentfunctionals can produce models with significantly different character, an essential step inthe solution of any inverse problem is to select the ‘right’ functional. Minimization of the 12norm of the model yields a ‘smallest’ model that often has high variability and hence maybe difficult to interpret. Similarly, models can be constructed by minimizing the 12 norm ofthe 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 asminimum structure solutions for many geophysical inverse problems. It is anticipated thatthe true earth model that gives rise to the measured data has at least as much structure asobserved on these constructed models.Methods for calculating minimum structure solutions have been described by variousauthors. These methods differ in the manner in which the required additional informationregarding the model is treated. Johnson and Gilbert (1972) construct a smoothest modelby explicitly incorporating endpoint values of the model and its derivative into the objectivefunctional to be minimized. Oldenburg (1984), in computing a flattest model, specifiesan a priori model value at an endpoint. Parker, Shure, and ilildebrand (1987) minimizea seminorm of the model. This chapter discusses two generalizations to the conventionalmethod for constructing a smoothest model. In this technique, the governing equations areintegrated by parts twice; each such integration introduces extra parameters into the modelconstruction problem. Numerical values for these additional parameters are obtained eitherby (i) prescribing a priori values at any abscissae within the interval of model definition,123or (ii) optimizing the appropriate objective function with respect to the parameters. Thesegeneralizations are readily extended to model construction problems that involve minimizingthe 12 norm of other derivatives of the model.In this chapter, the smoothest model construction formalism is applied to the problemof determining a depth proffle for a critically refracting horizon. It has only recently beendemonstrated that refraction traveltime data can be processed to yield point estimates ofinterface 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 methodassociated 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 locationalong the seismic line. The position of the refracting interface is then estimated by theenvelope of the set of arcs. This envelope is approximated by tangent lines to pairs of adjacentarcs; problems specific to this method are identified and discussed by Hatherly (1980). Thetechnique 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 relatedto calculating an envelope of a set of circular arcs. It also possesses several specific advantagesregarding treatment of the data. However, perhaps the most important difference pertainsdirectly to the derived models: the method proposed here yields the smoothest refractordepth profile consistent with the given data. This attribute is not shared by the modelconstructed via the GRM technique.A general theoretical developement of the smoothest model construction method is presented before discussing the particular application to interpolation. Hence, the results in thischapter are useful for the solution of other geophysical inverse problems that may require a124smooth model. Assume that an N-vector of observed data e0b3 is related to the true earthmodel mj(x) via a Fredholni integral equation of the first kind:be0b3= 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. Thetrue earth model mt(x) is unknown and is to be estimated. Predicted data generated by aconstructed earth model m(x) are given by the linear functionalbed(m) j g(x) m(x) dx. (6.2)In the following development, the 12 norm is used to quantitatively measure both modelstructure and data misfit. For example, the square of the 12 norm of the second derivativeof 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 isgiven byII e03 — e,,.d(m) 12 = [eobs —e1.d(m)j [cobs — elfl.d(m)},where the superscript T indicates transposition.6.2 Smoothest model constructionA detailed discussion of the method for calculating the flattest model is given by Aldridgeet al. (1991). The analogous derivation for the smoothest model is conceptually similar,although algebraically more complicated.1256.2.1 Modified data equationInitially, equation (6.2) is integrated by parts to obtain a modified data equation satisfiedby the first derivative of the constructed model:bed(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 ofequation (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 fundamentaltheorem of calculusb bm(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) yieldsed(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 ofthe abscissa Cl. This independence is emphasized by writing d1 for m(ci). Equation (6.6)becomesb,ed(m’) = d1h(b)—j p(x;cl)m’(x) dx, (6.7)126where a new kernel function vector is defined byp(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 derivativem’, 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 secondderivative of the constructed model can be derived. Thus, expression (6.7) is integrated byparts and the fundamental theorem is used to transfer endpoint information to anotherarbitrary location C2 with the interval [a, bJ. The result ised(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 byp(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) aredefined. 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 nonlinearfunctional of the second derivative of the constructed model. In the particular case whereboth c1 and C2 are restricted to the upper endpoint b, formulae (6.9) and (6.10) reduce tofamiliar forms.1276.2..2 Objective functionIn order to construct a smooth model that simultaneously minimizes the misfit betweenobserved and predicted data, consider the objective functionalmm” 112 + ‘N [eobi — eWd(m”)] 112 . (6.12)The scalar jt (0 < +oo) is a tradeoff parameter that controls the relative importance ofthe two terms. Equation (6.12) also includes a model derivative weighting function w(z) anda data weighting matrix W. The function w(z) can be used to emphasize certain portionsof 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 takento be the inverse of the covariance matrix of the observational error 6e. However, the matrixin (6.12) is not restricted to this specific type of weighting.Appendix D demonstrates that the second derivative function that extremizes 4 mustbe 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 straightforward to demonstrate thatwm” 112 = aTF(ci,c)a (6.14)ande,,.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):bF(c1,c2) j w(x)2p(x;cl,c2) p(x;ci,c2)T dx. (6.16)128The matrix F(c1,c2) is symmetric and, if the original kernel functions gj(x) are linearlyindependent 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 wellas explicit dependence on Cl). Note that expression (6.17) is a complete quadratic form inthe independent variables a, d1, and d2 in the sense that afl possible terms are present. Incontrast, 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 theconstant term. Additionally, the remaining terms differ in detail from those given above.6.2.3 Extremizing the objective functionStraightforward methods of multiva.riable calculus are now used to extrernize the objective function with respect to all of the variables relevant to the problem. Thus, calculatingO/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)129If values for the four scalars (Cl, d1, C2, (12) are prescribed, then (6.18) can be solved for thecoefficient vector a. This offers an interesting alternative to the conventional method forconstructing a smoothest model, where Cl and c2 are typically restricted to one of the endpoints a or b. However, for the current analysis, optimum values of the four parameters aredesired. 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 formsaTh(b)= 0, (6.19)aTk(b)= 0. (6.20)Geometricafly, these conditions imply that the N-dimensional coefficient vector a must beorthogonal 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 abscissaeCl 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 abscissaeCl 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, andd2, and can be solved simultaneously with standard techniques of linear algebra. However,130additional insight is obtained by first eliminating the coefficient vector a to form a simple2 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 aregiven byall = 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 CauchySchwartz inequality that the determinant of this system is nonzero. The optimum values ofd1 and d2 obtained by solving (6.21) are linear combinations of the observed data eobs. Afterthese two parameters are determined, the coefficient vector a is found by solving equation(6.18).6.2.4 Constructing the modelThe 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) reproducethe observed data e0b3 in the limit as the tradeoff parameter JL approaches zero. Thusm() = 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 beneficialimplications for the curvature of the model. In general, the second derivative of the model131is given by expression (6.13). Substituting in the explicit form for the kernel function vectorp(x;cl,c2) yieldsm”(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 tom”(z) = a Tk(z)/w(x)2 These extra conditions remove a step discontinuity in the secondderivative at a = cz due to the Heaviside function. Assuming that k(x) and w(z) arecontinuous, then the constructed model is twice continuously differentiable, i.e. m(x) E C2on (a,b).This model is also unique except in those cases where the determinant of the 2 x 2matrix A vanishes. Three particular cases are readily identified by analysis of det A. Fromequations (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 casearises in the following manner. Since () is positive definite, r—’/2(p) exists. Define twovectors u and v as:r’/2()k(b).The above expression for the determinant becomes det A = (uT v)2 — (uT u) (vT v). Ifu /3v (where /3 is a constant), then the Cauchy-Schwartz inequality implies that u 112v > (uTv)2. Hence, the determinant is nonvanishing in this situation, which isequivalent to the condition h(b) /3k(b).In these three cases, the constructed model is unique only to within an arbitrary additiveconstant ao or linear function a + aix. These situations are now described in detail.132Case 1: h(b) = 0. This situation arises if the original kernel functions g3(x) all possess zeroarea (e.g., Oldenburg, 1981). The system Ad = b reduces to a single equation that can besolved for d2:L-(ITp—1(d —e032 —— 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 thencalculated from (6.23) by picking any value for the parameter d1. Thus, constructed modelsdiffer by an arbitrary additive constant.Case 2: k(b)= 0. The 2 x 2 system Ad b reduces to a single relation between d1 andd2: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 aboverelation between d1 and d2. Constructed models then differ by an arbitrary linear functionof x.Case 3: k(b) = /3 h(b). This case is a generalization of the previous one. The two equationsin (6.21) are no longer independent, but reduce to the single expressiond1 + (b—ci —/3)d2 =133Once 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 2above. Hence, the model is unique only to within an arbitrary linear function.This analysis highlights a particular advantage of the current formulation, where theseparate system (6.21) is derived for the two parameters d1 and d2. The specific mathematical conditions for nonuniqueness in the calculated model can be determined simply byexamining the equation det A = 0. Moreover, as indicated above, model construction canstifi 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 modelThe theory developed in the previous section is clearly illustrated by the problem ofconstructing a smooth interpolation of a set of discrete samples. This is a common problemin many branches of geophysics. The particular application considered here consists ofcalculating a continuous depth profile of a refracting interface from a set of point estimatesof the refractor depth and dip. These point estimates can be derived by analyzing thefirst break traveltimes of a seismic experiment by various techniques. There are many suchtechniques available to the practicing interpreter, and these are not examined individuallyhere. However, the methods typically employ the assumption that the critically refractinghorizon is plane or nearly plane. Under these circumstances, it is logical to construct aninterface depth proffle that possesses minimum curvature. Thus, the final model for theinterface adheres closely to the prior theoretical assumptions used for inferring its depth.Many shallow seismic refraction projects are undertaken in conjunction with a driffingprogram (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 refraction134method is then used to extend this information between or beyond the wells. In somecases, outcrops of the target horizon provide additional geologic control on both depthand dip (Kilty et at, 1986). When constructing a proffle of the interface, the depth anddip information from all available sources should be integrated together. Extra data fromdri]lholes, trenching, outcrops, etc. then contribute directly to the solution, rather thanmerely providing ‘tiepoints’ to the profile determined from refraction data. This goal can bereadily achieved if the problem of caiculating the refractor proffle is posed as an interpolationissue. In this case, the data weighting matrix W assumes an important role. This matrixcan be used to nondimensionaJize different data types, and to emphasize those particulardata deemed more important in the inversion. In contrast, there is no provision in theGRM depth conversion technique for including additional geological or geophysical datawith variable weighting.6.3.1 General theoryThe complete set of N measured data is divided into two distinct subsets: n depthsamples and N — n dip samples, where 0 m N. These data are observed at abscissae ajthat are ordered as follows:Depth: a < X2 < < < xDip: a Xn+1 < Xn+2 < < XN1 < XN b.The abscissae associated with the dip samples do not necessarily coincide with those of thedepth 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 aretreated first. For j = 1,2,.. . , n, a predicted datum is considered to be a sample from a135function e(x):eTd e(xj)=The th component of the kernel function vector g(x) is a Dirac delta function centered atthe sampling point 2j : gj(z) = S(x — z,). Hence, the singly and doubly integrated kernelsareh(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 redefinedas slope, rather than dip samples. The transformation is easily effected by calculating thetangent of each dip angle. From the fundamental theorem of calculus, the difference in theslope of the function e(x) between two locations z1 and C2 is given byZj &e’(xj)—e’(c2)= £2e”(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 secondderivative of e(x):e’(x) = d2 + jbpi(z;cl,c2)eII(x) dx,where d2 = e’(c2) and the kernel function is defined bypj(x;cl,c2) = H(x—c)—H(x-—xj). (6.25)136Comparison of the above two expressions with equations (6.9) and (6.10) indicates thath(b) = 0 and k(b) = —1 for the index range j = n + 1, m + 2,. . . , N. Hence, if the datasetconsists entirely of slope samples (n = 0), then the interpolation problem is an example ofthe Case 1 phenomenon described in the previous section (i.e., the constructed interpolante(x) is unique only to within an additive constant).If the weighting function w(x) equals unity, then explicit formulae for the elements ofthe inner product matrix F(ci, c2) can be derived by substituting the above kernel functionsinto (6.16) and integrating. Arbitrary nonuniform weighting requires that the integrals beevaluated 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 tovarious locations of the coordinates z and x relative to c1 and C2. These equations areomitted 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 beused.As indicated previously, it is possible to solve for a, d1, and d2 simultaneously. This isthe preferred method of solution for situations with a large amount of data. It eliminatesthe need to invert F(t) in order to obtain the elements of the matrix A and the vector bas in (6.22). However, for small scale problems like the following examples, calculation ofF(t) is not computationally burdensome. Singular value decomposition (SVD) is usedhere to compute the inverse. An advantageous feature of this technique is that it allows thesingular vectors associated with very small eigenvalues to be winnowed from the solutionsimply by truncating the SVD sum. Small eigenvalues arise when the matrix F(Au) is nearlysingular, 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+l137These quantities are all defined for the full index range i = 1,2,. . . , N. After some algebraicmanipulation, the elements of the matrix A and the vector b in equation (6.21) reduce tothe simple formsall P, a12 = (xj — ci )P +i=l 2=1 i=n+ln n Na21 (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 completed by constructing the function e(z) via equation (6.23). Expressions for the doublyintegrated 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 todifferent location of the coordinate x relative to the two abscissae Cl and C2.6.3.2 Numerical exampleFigures 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/bedrockinterface (dashed line). A small channel is located adjacent to a larger and broader anticine.Elevations are referenced to an arbitrary horizontal datum plane. Nonuniformly spacedelevation and slope samples from the interface are indicated by asterisks. Initially, these areconsidered to be error free.If accurate estimates of interface depth and dip are available (say, from a shallow borehole) then these values can be used to construct a smoothest elevation model. Figure 6.1depicts this situation. The prescribed elevation and slope values are indicated by open138—020 10 20 40 50‘ 0.20 :50x(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 andslope are taken from the true earth model and are indicated by small squares at abscissaci = C2 = 23.75 m. The 12 norm of the second derivative of the constructed model is0.228 m_h/2.6.00 100240 501390.20.0—020 10 20 0 40 500.10_v—o.io0 10 20 30 40 50x(m)Fig. 6.2. Same as Fig. 6.1 except that optimum values of the refractor elevation and slopeare calculated at Cl = c2 = 23.75 m. The 12 norm of the second derivative is reduced toII e” = 0.129 m’2.6.00 10 40 50-II I I140squares at horizontal coordinate ci = C2 = 23.75 m. The constructed model and its derivative (solid curves) pass through these points. However, the second derivative of the modelpossesses a jump discontinuity at this location. This discontinuity arises from the Heavisidestep 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 thediscontinuity in e”(x) can be magnified (Aidridge et al., 1991).Figure 6.2 displays the smoothest elevation model constructed by using the optimumvalues of d1 and d2. There are obvious differences from the previous model on the flank ofthe anticine and on the flat portion between anticine and syncine. The second derivative ofthe elevation proffle is now a continuous, piecewise linear function. It is emphasized that thesame model is obtained if the parameters d1 and d2 are optimized at any pair of horizontalpositions ci and C2 along the proffle (including outside of the range encompassed by thesix elevation samples). The 12 norm of the second derivative is reduced by optimizing theparameters, as expected. Aidridge et at (1991) discuss the relation of these interpolants tocubic splines. In particular, the constructed elevation profile in Figure 6.2 is a cubic splineon [a, bJ, whereas e(z) in Figure 6.1 consists of two cubic splines joined end-to-end at thecoordinate x = C2.Inclusion of slope samples, as in Figure 6.3, results in a more accurate reconstruction ofthe 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 derivativeat each abscissa x associated with a slope datum. Hence, the constructed elevation proffleis no longer a cubic spline on [a, bj.In actual practice, the calculated interface depth and dip samples will contain somerandom error. This situation is simulated in Figure 6.4. Random numbers drawn from auniform probability distribution on ±0.50 m (standard deviation = 0.29 m) are added tothe accurate elevations. Also, each slope sample is perturbed by an amount correspondingto a uniformly distributed random dip angle on ±5° (standard deviation = 2.9°). A modelthat is an exact fit to these error contaminated samples would contain spurious structureinduced by the noise. Hence, a model is constructed that reproduces the erroneous data1410.20‘.0.10V—0.1050Fig. 6.3. A smoothest refractor model constructed from six depth samples and six dipsamples. All data are error-free. Optimum values of the refractor elevation and slope arecalculated at Cl = C2 = 23.75 m. e” 0.354 m’2.6.0‘p0 10 20—0.2400 10 20 40 500 10 20 30 40x(m)1427.016.0500.2nN-‘-‘ 004,—02-4.0.10N -____v—0.100 ‘10 20 30 40 50x(m)Fig. 6.4. A smoothest refractor model constructed from error contaminated data. Depthand dip samples are located at the same abscissae as in Fig. 6.3. Optimum values of therefractor elevation and slope are calculated at c = C2 = 23.75 m. J= 0.259 m’2.0 tO 20 50 400 tO 20 SO 40 50—1’-•1 ‘143approximately, rather than exactly. The tradeoff parameter ,u and the weighting matrix Wcan be used advantageously for this purpose. Since a constructed smoothest model dependson the value of the tradeoff parameter, predicted data generated by this model will also bea function of p. The notation ez,,.d(p) is used to explicitly denote this dependence. If theweighting matrix is diagonal with elements W:, i = 1,2,. . . , N, then the degree of misfitbetween observed and predicted data is quantified byEdepth(JL) j,2 {e — erd(p)]2,6slope(p) N —W [ebs — e()j2.- jfl+1A value for the tradeoff parameter is sought such that these misfits are approximated bydepth(P) Eslope(P) N—nNote that, for uniform weighting, the right hand side of each expression above is proportionalto the rms value of the standard deviations of the data errors. Selection of a search procedurefor p is largely a matter of personal preference. For small scale problems like the presentexample, a few trial runs with guessed values of p allow a sufficiently precise estimate to bemade. Alternately, it is possible to implement an automated iterative scheme that convergesupon 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 uniformdistribution of dip angles on ±5°. Also, optimum values for the two parameters d1 andd2 are used. Interest in detecting the shallow channel may stem from various exploration144objectives: placer ore localization, groundwater accumulation, contaminant migration, etc.However, the figure suggests that the channel is just at the limit of detectability for thegiven horizontal sampling and error statistics.6.4 ConclusionThe conventional method for constructing the smoothest model is generalized in twoimportant ways:1) The additional information required to calculate the smoothest model may be specifiedanywhere within the interval of definition [a, bj. This improvement can be particularly usefulif 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 andits slope directly frorri the observcd. data. This technique is advisable if independent valuesof these parameters are unknown or are difficult to estimate accurately.It is evident that these techniques can be extended to model construction situationswhere the 12 norm of a higher derivative of the model is minimized. The modified dataequation satisfied by the m derivative of the model m()(z) will entail ii parameters cij defined at the abscissae cj as follows: dj m(1)(cj), i = 1,2,... ,n. Extremizing the relevantobjective function leads to an n x n system of equations Ad = b for the unknown d. Aunique solution exists when det A 0, and the model constructed with these parameterswill possess a continuous th derivative (m(a) E C’s) provided that the original kernel function 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 towithin an additive polynomial of at most degree n — 1. Only the smallest model (n = 0) isunique in all situations.Application of this model construction formalism to the problem of calculating an interface depth profile yields a useful tool for seismic refraction exploration. The smoothestmodel is the natural model to adopt in this situation because refraction traveltime inversion145methods assume, either explicitly or implicitly, that local interface curvature is negligable.The resulting depth profile and its derivatives are mathematically well defined throughoutthe horizontal range (a, b) and thus can be used for other forward modeling purposes. Themethod is also flexible, as evidenced by the following specific advantages compared withdepth determination via the GRM:1) Additional depth and clip data arising from a variety of geological, geophysical, or engineering 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 feasibility of the technique for refractor depth profile construction. Obviously, there is stifi muchto learn about issues such as weighting of different data types, resolvability of small scalefeatures, optimum horizontal sampling, etc. However, the algorithm presented here providesa flexible tool for investigating these and related phenomena.146CHAPTER 7REFRACTOR IMAGING USING AN AUTOMATEDWAVEFRONT RECONSTRUCTION METHOD7.1 IntroductionThe Wavefront Method is one of the earliest of the many techniques for interpreting refraction arrival times. In 1930, Thornburgh demonstrated that subsurface wavefronts couldbe reconstructed from surface arrival times by applying Huygens’ principle in reverse. Subsequently, Hagedoorn (1959) elucidated an imaging condition for delineating a refractinghorizon. First, two oppositely propagating wavefront systems are reconstructed from thearrival times recorded on a forward and reverse spread, respectively. Then, pairs of thesesubsurface wavefronts intersect on or slightly below the refracting interface when the sumof their times equal.s the known reciprocal time (the shot-to-shot traveltime). This imaging principle yields the correct spatial locus of a critically refracting horizon if the earthconsists of constant velocity layers ‘bounded by plane dipping interfaces. However, severalinvestigators have demonstrated that the imaging condition is reasonably accurate even ifthe 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 fromnonpiane structure on the refracting interface, or a velocity gradient within the underlyingmedium.Extensive application of the wavefront method has been limited by two factors: i) laborious 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 addresses the first of these two issues. Instead of defining the wavefronts by a tedious graphicalapplication of Huygens’ principle (e.g., Rockwell, 1967), a finite-difference computer algorithm 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 velocitystructure.147Recently, Hill (1987) downward continued refracted waveforms to obtain a two dimensional image of shallow structure. Although the technique presented here works only witharrival times, the goal is identical. The advantage of this approach resides in its computational simplicity. Since the propagation algorithm operates directly in the space-timedomain, no transformations of the recorded wavefield, with attendent concerns about sampling adequacy (Clayton and McMechan, 1981; Hill, 1987) are necessary. Furthermore, trueamplitude recording and processing of seismic traces are not required. However, prior picking of these traces to obtain the arrival times is necessary, and this may be a time consumingjob in some situations.7.2 Finite-difference traveltimes7.2.1 Wavefront ConstructionVidale (1988, 1990) has recently developed an algorithm for calculating the first arrivaltimes 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-differenceoperators 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 algorithm properly handles the various wave types that comprise first arrivals (body waves,head waves, and diffractions). Subsequent contouring of the computed traveltime field yieldsa visual impression of propagating wavefronts. The phrase wavefront construction is usedhere to refer to traveltime loci calculated in this manner. Figure 7.1 depicts the subsurfacewavefront systems generated by a sequence of shots buried in an earth model with undulating surface and refractor topography. The direct wave through the overburden is the initialarrival near each shot location. Beyond the crossover distance, the wave refracted by thehigher velocity bedrock arrives first. The traveltimes recorded along the nonplane surface ofthe model are accurately computed by assigning a P-wave velocity to the uppermost layerequal to the speed of sound in air (“-i 350 m/s). Wavefront contours are then suppressed indepth(m)2010depth(m)0302010300 a’0300depth(m)20100300depth(m)depth(m)201003020100—CnC12e4-.,—.b.•c-,-.-CCDII.C)’e4I)÷:z•QCD‘-CD )-jp, I-,oC•-D-0 0CDCDp,)-CDCD.<CI)CDqCDC’)Cn_I-+-.CD,0C.)$CDCDCDCDIII-’.CI:,Cl) Eo—--0 N 0a’(I) I-.0La (7’ 0’ 0 -1 0 I 0 0I 0 0-2 011490t/00horizontal 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 illustratethat nonpiane topography has a complicating effect on an interpretation.Vidale’s wavefront construction algorithm has been altered in two important ways inorder to improve its suitability for the shallow refraction problem. First, the traveltimecalculations are initiated from a spatially extended source, rather than a point source. Inthe 2D case, source activation times are specified on and inside a rectangular region locatedwithin the velocity field; arrival times at grid points outside this rectangle are generated bythe normal working of the algorithm. Since wavefronts are strongly curved in the immediatevicinity, of a point source, use of the plane wave finite-difference operators will yield inaccurate traveltimes in this region. Moreover, these inaccuracies will be propagated to greaterdistances, where the plane wave extrapolators are locally valid. In order to avoid this problem, near source traveltimes are calculated via mathematically exact formulae appropriate0 25 50 75 100150for either a constant or linear velocity field. Although more complicated velocity distributions can be considered, these particular velocity functions provide sufficient flexibility formany traveltime computation problems.Second, the mathematical form of the traveltime extrapolation operator is modified inthose cases where there is a large velocity increase across a grid cell. This situation is relatively common in the shallow refraction environment. The interface between unconsolidatedoverburden and consolidated bedrock, or between saturated and unsaturated alluvium, oftenrepresents a sharp velocity increase. In these cases, as the following analysis indicates, theconventional traveltime extrapolation formula may fail.Figure 7.3a depicts a system of plane wavefronts propagating across a square grid cellwith side length h. The arrival time at the corner numbered 4 must be calculated from theknown arrival times t1, t2, and t3 at the other three corners of the cell. Assuming a planewave advancing with a constant slowness s, this time is given by t4 — ti + (/h coswhere the angle S describes the ray direction relative to the cell diagonal. Simple geometricanalysis yieldsI t3_t2COSS = 1—____Hencet4 = t1 + i/2(hs)2 — (t3 —t2) . (7.1)This expression is identical to Vidale’s (1988) equation (3), which was derived by approximating the partial derivatives in the 2D eikonal equation by finite differences and thensolving algebraically for t4. The present derivation clearly reveals the underlying geometricassumption of plane wave propagation.If the argument of the square root in equation (7.1) becomes negative, then the planewave extrapolation formula is obviously invalid. This may occur, for example, if there is adramatic velocity increase across the cell (implying that the slowness s assigned to the cell151Fig. 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 planewavefront approximation is abandoned and Huygens’ principle is used directly to calculatethe next traveltime. Although this computed time is not exactly correct, extensive numericaltesting 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 t2b)ti t2ht3 t4t3 t4\\1527.2.2 Wavefront reconstructionFigure 7.4 displays the forward and reverse wavefront systems generated by shootingover a shallow syncine. These are the wavefronts that give rise to the first arrival timesobserved on the surface. The finite-difference traveltime algorithm can now be used torecreate 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 horizontaldimension. This line source is then activated sequentially (rather than simultaneously) withan 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 recordedfrom distant shotpoints (Rockwell, 1967; Ackermaun et al., 1986). The line source generatesa set of wavefronts radiating downward into the specified velocity field (Figure 7.5). Thedownward continuation velocity function v(z, z) is selected as a good approximation tothe actual near surface velocity structure. Hence, within the overburden, the calculatedwavefronts coincide with the emerging refracted wavefronts of Figure 7.4. Since the positionof the refracting interface is initially unknown, the wavefronts are continued to greater depthusing the known velocity field v(c, z). Rockwell (1967) referred to these traveltime loci asa “directed wavefront system”; in this study, the phrase wavefront reconstruction is used todescribe the process of creating an emergent wavefront system from recorded surface arrivaltimes.7.3 Refractor imagingLet tf(x, z) and tr(z, z) refer to the subsurface traveltime fields reconstructed from theforward and reverse arrival times, respectively. Then, according to Hagedoorn’s (1959)imaging principle, the refracting interface is implicitly defined by the relation1532.5Ici)2.5Fig. 7.4. Forward and reverse first arrival wavefronts for a shallow syndine model. Overburden velocity v = 1500 m/s, bedrock velocity v2 = 2500 rn/s. Grid cell size is 5 m andcontour interval for wavefronts is 30 ms.0.501 1.5 2.5 1 1.5horizontal position (km)2154Ia)2.5Ia)U)•0Fig. 7.5. Emergent wavefronts reconstructed from the surface arrival times recorded overthe shallow syndine. Downward continuation velocity v(x, z) = 1500 rn/s. Grid cell size is5 m and wavefront contour interval is 30 ms.0.501 1.5 2.5 1 1.5 2 2.5horizontal position (km)155tf(X,Z) + tr(x,z) TR. (7.4)Figure 7.6 graphically illustrates the superposition of the two reconstructed wavefront systems shown in the prior figure. This 2D array of superposed traveltimes is systematicallysearched 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 usedto find the proper depth. The resulting depth locus z() (dashed line in Figure 7.6) is anaccurate spatial image of the original refracting horizon, except near the edges of the inputvelocity field where the subsurface wavefronts are not reconstructed correctly.Figure 7.7 indicates that the technique is also capable of imaging anticinal structureThe apex of the anticine is imaged slightly too deep because the refracted rays penetratebeneatk 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 ofsyncinal structure.The calculated locus for the refracting horizon depends on the reciprocal time TR andvelocity field v(x, z) used for downward continuation of the surface arrival times. Variationsin these quantitites from their correct values will induce variations in the depth and positionof the refractor.It is relatively easy to assess the dependence of the refractor image on the value of thereciprocal time. The forward and reverse subsurface traveltime fields are added togetherand the result is contoured for various candidate “imaging times”. Figure 7.8a illustratesthis situation for the buried syndline. If the imaging time used is less than or greater thanthe 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 bythe classical wavefront method (Rockwell, 1967, p. 378). The difference arises from themethod of reconstructing the subsurface wavefronts. Currently, finite-difference traveltimecomputations are initiated with the source function (7.3) and the algorithm is run forward1560IU)0 2.50I1)0 .5 1 1.5 2 2.5horizontal position (km)Fig. 7.6. (a) Superposition of the two reconstructed wavefront systems of Fig. 7.5. Dashedline is the locus satisfying the refractor imaging condition. (b) Comparsion of the true (solid)and imaged (dashed) refracting interfaces.0 .5 1 1.5 2I I I I b)I I1570I2.50I.4)U,C’20 .5 1 1.5 2 2.5horizontal position (1cm)Fig. 7.7. (a) Superposition of the forward and reverse reconstructed wavefronts (contourinterval = 30 ms) over a shaJiow anticine. Dashed line is the refractor image. (b) Comparisonof the true (solid) and imaged (dashed) refracting interfaces.U).5 1 1.5 2I I I Ib)I I I1580I__________________Lt)C2.5 1 1.5 2 2.50Ir-Li)e2.5horizontal position (km)Fig. 7.8. Dependence of the syncline locus on imaging time (top) and downward continuationvelocity (bottom). Different images correspond to increments of 10 ms in imaging time and100 rn/s in velocity, respectively. Heavy lines are the images corresponding to the correctvalues of TR and v.I I I___ __ ____TR—6O ms1TR+6O ms0 .5 1 1.5 2159in time. Hence, subsurface wavefronts are labeled with times later than the surface sourcevalues. In contrast, the classical wavefront reconstruction methods label the subsurfacewavefronts with times earlier than the surface measured times. In either case, the trueposition of the interface corresponds to an imaging time equal to TR.Quantifying the dependence of the refractor image on the downward continuation velocity is more complicated. Forward and reverse subsurface wavefront systems must bereconstructed for each velocity function used in the analysis. Figure 7.8b displays a set ofimages of the shallow syncine calculated for various values of a constant downward continuation velocity. If the velocity is less than or greater than the actual overburden velocity, thenthe interface image is too shallow or too deep, respectively. Moreover, a grossly incorrectcontinuation velocity distorts the shape of the interface structure. Hence, in common withmany other seismic refraction interpretation techniques, accurate time-to-depth conversionwith the wavefront method requires good knowledge of the overburden velocity distribution.This information can be obtained from uphole times, direct and reflected arrivals, shallowrefractions, and borehole data.Finally, the accuracy of the solution depends on the reliability of the picked first arrivaltimes. The ability of the method to resolve small scale features on the refracting horizon isalso limited by the field geophone interval. These phenomena are analyzed by performingthe inversion with noisy traveltime data sampled at an assumed geophone interval. Figure7.9 displays the refractor image obtained by downward continuing error contaminated arrivaltimes sampled every 25 m. Spatially correlated, normally distributed time errors (standarddeviation 5 msec; correlation distance = 100 m) are added to the theoretically exactrefraction picks. Appendix F describes the method used for computing correlated randomnoise. 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 thesyndine, and its lateral position and depth are approximately correct. Long wavelengthundulations on the refractor are artifacts of the spatially correlated noise, but are not undulyharmful to the structural interpretation.1600I4-’LI)0 2.5Fig. 7.9. Reconstructed wavefonts and syndine image formed from noisy traveltime datasampled at a geophone interval of 25 m. Grid cell size is 5 m and wavefront contour intervalis 30 ms. Long wavelength undulations on the refractor arise from spatiafly correlatedtraveltime errors.7.4 Refractor velocity estimationA particular advantage of the wavefront method is that the interface depth calculationis independent of the refractor velocity. Rather, the velocity of the substratum can be estimated after the position of the refracting horizon is determined. The distance betweentwo points on the interface divided by the difference in the reconstructed wavefront times atthese points is an estimate of the refractor velocity. This value is assigned to the midpointof the two points for plotting purposes. In effect, the directional derivative of the subsurfacetraveltime 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 forward or reverse wavefront systems may be used for the computation. Figure 7.10 illustratesthat the refractor velocity estimated in this manner possesses systematic errors related to.5 1 1.5 2horizontal position (km)161CO(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, shortwavelength oscillations are artifacts of the grid interval.Coc’i.75 1 1.25 1.5 1.75.75 1 1.25 1.5 1.75horizontal position (km)162the interface structure. However, for the two synthetic examples examined here, the inferredvelocity values are everywhere within 3% (or ± 75 m/s) of the correct values.7.5 Field data exampleThe interface imaging procedure is tested with a shallow refraction dataset acquiredat the archeological site of Phalasarna in western Crete. Hadjidaki (1988) discusses thehistorical and archeological significance of this site and also gives a detailed description ofthe surface and near subsurface conditions. Forward and reverse refraction proffles wererecorded along an inline spread of 18 geophones (geophone interval = 0.5 m, near sourceoffset 0.5 m) during the summer of 1989. The data acquisition system consisted of aportable signal stacking seismograph with a hammer energy source. First arrival time pickswere 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 spreaddiffer 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, reciprocaltime corrections (Hatherly, 1982) are applied to the picked arrival times. A constant timeshift is added to the raw time picks on each source gather in order to adjust the observedreciprocal times to the average value of 16.55 msec.Figure 7.11 displays the 18 first break picks recorded on the forward and reverse profflesafter application of these reciprocal time corrections. A preliminary interpretation of theplotted 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 surfacevelocity function for subsequent downward continuation of the refracted arrival times. Acubic spline is fitted to the 16 refraction picks on each spread and is extrapolated to zerooffset as a straight line. These curves are then used in equation (7.3) to calculate the sourceinitiation functions required for the wavefront reconstruction algorithm.163horizontal 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 bysmooth curves.Figure 7.12a depicts the subsurface wavefronts generated by downward continuing therefracted 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 correctedreciprocal time TR = 16.55 msec. The refractor velocity estimate (Figure 7.12b) exhibitstwo distinct zones: i) abrupt variations about 800 rn/s on the left, and ii) a low velocityzone slower than 800 rn/s on the right.The validity of various refractor velocity functions can be tested by using the wavefrontconstruction algorithm to compare predicted traveltimes with the observed traveltimes. Auniform 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 andreverse wavefronts of Figure 7.14 between 6 m and 8 m horizontal position. Note that00 1 2 3 4 5 6 7 8 9r -.-.Ct,CD I-.0 ‘-4‘-i.--CD.CDCDç”aqU)CD0..I4.9CD--‘%—,o CDa.,. 0ojc4•4-”O 0—‘01-sU) 0 I-.o CD 0U)-I.velocity(m/s)200400800800100002depth(m)100CA)Ci)0 —.N 0 0 (1) c-I 0Q1Q1-2-2coCl)-4-)C)0G.)CD‘-4a)0CD‘-IVi‘-4F000‘-I00CD00CD0000CQ01651 2 3 4 5 6 7 8 90 1 2 3 4 5 6 7 8 90horizontal position (m)Fig. 7.13. Comparison of observed and predicted first break times for two refractor velocityfunctions. (a) Constant velocity (dashed) equals 800 rn/s and variable velocity (solid) includes a low velocity zone between 6 and 8 m. (b) Predicted traveltimes calculated with theconstant refractor velocity. (c) Predicted traveltimes calculated with the variable refractorvelocity..4-,166Fig. 7.14. Subsurface wavefronts (contour interval = 0.5 ms) constructed from the earthmodel with the laterally varying refractor velocity. Grid cell size is 0.02 m. Note the effectof the low velocity zone between 6 and 8 m.0.4-,C)0 102 3 4 5 6 7 8 90 1 2 3 4 5 6 7 8 9horizontal position (m)167the strong variations in refractor velocity displayed on the left in Figure 7.12b need not beincorporated into the model to obtain an adequate fit to the measured arrival times. Furtheradjustment 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 interpretationsof the observed traveltime data, incorporating multiple layers or lateral changes in structureand/or velocity, are possible. Since the refraction dataset does not include arrival timesrecorded 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 subsurfaceinformation from the vicinity of the refraction proffle. Archeological treadling conductedabout 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, aeratedsandstone might have a P-wave velocity as low as “.‘ 800 rn/s. Hence, the preliminaryinterpretation is that the interface imaged in Figure 7.14 is the upper surface of the sandstonebedrock. The very low velocity between 6 m and 8 m may be a zone of more extensiveweathering, fracturing, or aeration. Alternately, it is possible that one of the overlyingshallow gravel layers has been imaged.7.6 ConclusionThe 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 rapidcalculation of the subsurface wavefront systems that give rise to the recorded arrival times.168Although the synthetic examples described here use a uniform near surface velocity, downward continuation through a varying velocity field is also possible with no increase in computation time. The buried refracting horizon is delineated in a subsequent step by applyingHagedoorn’s imaging principle. No prior assumption regarding the refractor velocity is required. Rather, the velocity of the substratum can be estimated by calculating the directionalderivative of the reconstructed wavefront systems along the imaged interface.Picking of first arrival times and assignment of these picks to specific refractors arenecessary 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 velocityfield. However, since the technique is not computationafly intensive, it is possible to assessthe magnitude of the position and depth uncertainty by performing the inversion repetitively. The forward modeling capabffities of the finite-difference traveltime algorithm canalso 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 thesignificance of various features of the recovered earth model.Finally, two specific problem areas with the automated wavefront reconstruction methodhave been identified that merit further research: i) downward continuation of traveltimesrecorded along a nonplane surface, and ii) correction of the refractor velocity function forthe effects of structurally induced errors. Although a fully automated solution to theseproblems is not yet available, this should not prevent the immediate application of themethod in shallow seismic refraction exploration.169CHAPTER 8TWO DIMENSIONAL TOMOGRAPHIC INVERSION WITHFINITE-DIFFERENCE TRAVELTIMES8.1 IntroductionCurved ray traveltime tomography was originally developed by Bois et al. (1972) for thepurpose of estimating the seismic velocity distribution between two boreholes. Followingtheir seminal work, several investigators advanced the technology of tomographic imagingwith 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 accuratelyreconstructing the velocity field in highly refractive media. Although straight ray techniquesare adequate in media with relatively small velocity variations, it is difficult to decide whenthe simplifying assumption of straight line raypaths becomes invalid. Hence, there is acompeffing reason to use curved ray methods in all situations: they are based on a moreaccurate model of wave propagation through variable velocity media.Traveltime tomography is a nonlinear inverse problem that can be solved by local linearization and iteration. Since the velocity field is updated on each cycle of the tomographicimaging procedure, rays have to be retraced between all source-receiver pairs. This raytracing 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 mayencounter difficulties due to shadow zones and multipathing. In this study, the problemsassociated with conventional raytracing are circumvented by calculating traveltimes to allpoints of a two dimensional slowness field with a rapid finite-difference algorithm (Vidale,1988). Raypaths are then generated by following the steepest descent direction throughthe computed traveltimes from each receiver back to the source. This method yields theraypaths of all wave types that comprise first arrivals (body waves, head waves, and diffractions). Moreover, since arrival times are calculated throughout the slowness field, arbitraryrecording geometries are easily accommodated.170In addition to rapid and accurate forward modeling, tomographic inversion requires thesolution of a system of linear algebraic equations to obtain the improved velocity field. Thissolution may exhibit erratic and unphysical behavior due to noise in the observed traveltimedata and/or ill-conditioning in the equations. Hence, some form of regularization is oftenused to stabilize the inversion. For example, Lytle and Dines (1980) introduce Laplaciansmoothing into the system when calculating a perturbation to a slowness model. Bishop etal. (1985), Bregman et al. (1989), and White (1989) limit the size of the model perturbation by using the damped least squares method. Macrides et al. (1988) impose inequalityconstraints on a perturbation calculated via an ART algorithm. The present approach isto apply linear equality constraints directly to the slowness model, rather than to a modelperturbation, on each iteration of the inversion procedure. In addition to improving themathematical conditioning of the system, the constraint equations allow the introduction ofa priori geological or geophysical knowledge about the model into the iuversioa. In particular, the constraints may arise from a desire to impose a preferred character, like flatness orsmoothness, on the slowness solution. Alternately, one may seek a model that is close, insome quantitative sense, to a prescribed base model. Inclusion of these constraint equationsrestricts the nonuniqueness that is common in realistic tomographic inverse problems.The nonlinear tomographic inversion procedure described in this chapter consists of fourbasic steps:1) calculation of first arrival traveltimes from each source location to all points of a griddedslowness 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 existingslowness model,4) updating and (optionally) smoothing the slowness model.This four-step process is initiated with an estimate of the true slowness function, and isrepeated until an acceptable match is obtained between observed and calculated traveltimes.171The initial estimate is usually a uniform slowness field. Subsequent sections describe thesesteps in more detail, and demonstrate the inversion procedure using synthetic traveltimedata from simulated VSP and crosswell experiments.8.2 General theory and model representationThe traveltime of a seismic wave propagating through a slowness field s is given by thepath integralt(s)= f sdl, (8.1)r(s)where r(s) is the raypath connecting source and receiver, and dl is an incremental pathlength. Since the raypath locus depends on the slowness, the traveltiine t(s) is a nonlinearfunctional of s. The problem can be linearized by considering the traveltime differenceAt t(s + As) — t(s), where As is a perturbation to the slowness field s. For a sufficientlysmall 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 fieldcan be evaluated by integration along the unperturbed raypath. The traveltime differencebecomesAt= 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 raypathis a straight line segment. For a collection of raypaths, equation (8.2) is expressed as thematrix/vector productAt = A(m) Am, (8.3)where A(m) is a matrix of raypath length segments within the square cells of the slownessmodel m. In the tomographic inversion problem, a model perturbation Am is sought such172that the improved model m + m approximates the true slowness model that generates themeasured traveltime data. Hence, the traveltime difference vector it is given byt0b3 t,.d(m), (8.4)where t3 is a vector of observed arrival times, and t,.d(m) is a vector of predicted trayeltimes computed from the known slowness model m. In principle, equation (8.3) can besolved for the required model perturbation. In practice, solution difficulties arise because the-raypath matrix is commonly nonsquare, large, sparse, and rank deficient. Moreover, sincethe observed traveltimes contain random errors, the system (8.3) may be inconsistent. Inthis 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 iterationsof the model updating procedure may be necessary before the inagnitude of the traveltimeresidual vector zt becomes acceptably small.8.3 Forward modelingThe first arrival times of a seismic wave propagating through a two dimensional velocity structure are computed by Vidale’s (1988) finite-difference scheme. This algorithmuses plane wavefront traveltime operators to extrapolate arrival times from point to pointthroughout a uniformly spaced grid. The method is rapid and accurate, and can be applied to a heterogeneous medium with moderate to strong velocity variations. Podvin andLecomte (1991) describe improvements to the local traveltime extrapolators that allow models with very strong velocity contrasts to be examined. The traveltimes of all wave types thatcomprise first arrivals (body waves, head waves, and diffractions) are calculated. Reflectionsand other later arrivals are not included; this represents a limitation of the technique ascurrently formulated.Vidale’s method is based on a centered finite-difference solution of the eikonal equationon each square cell of a gridded slowness field. Thus, the associated discretization erroris second order in the grid cell size. An input slowness function s(x, z) is sampled on a173uniformly spaced two dimensional grid. If the grid interval is h, then the sampled slownessvalues are given by sj = sQc, zj), where ij = Xmin + (i — 1)h and Z = Zmjn + (j 1)h (withi 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 slownessvalues at the four bounding grid points:1mk = (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 vectorconstitute a row-ordered sequence of the two dimensional array of cell slowness values. Anindividual cell is referenced either by the coordinate indices of its upper left corner (ii), orby its sequential index (k).Calculations are initiated at a designated source point (,z3) (not necessarily coincident with a grid point) within the slowness model. Since wavefronts are strongly curvedin the immediate vicinity of a point source, plane wavefront traveltime extrapolators areinappropriate in this region. Furthermore, near source inaccuracies are propagated to allgreater distances. In order to mitigate these effects, the traveltimes in a near source rectangle are calculated via mathematically exact formulae derived from certain simple velocitydistributions. Hence, consider the linear velocity functionv(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 verticalcomponents of the velocity gradient vector. Numerical values for these constants are obtainedby 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,, + aand tan 4 = a,,/a. If the frame of reference is rotated through the angle , then the velocity174field (8.6) is transformed into a one dimensional function. The spatial coordinates in therotated system are designated by primes and are given byxl = 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 arrivaltime at an arbitrary near source location (z’, z’) can then be calculated by standard 1Dtechniques. The result isI 21( i l’2j( I\2I l — — a LkX Xs) IZ Z31LIx ,Z — Siflit ii ,a y 4v[v3+a(z —z3)JWavefronts associated with this time field are eccentric circles with centers that are displacedalong 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 theproper expression for a uniform velocity field (i.e., circular concentric wavefronts).The finite-difference algorithm calculates a traveltime tjj at every grid point of theslowness field. If a receiver is not located on a grid node, then an interpolator is neededto estimate the arrival time at the actual receiver position. Simple bilinear interpolationprovides adequate accuracy. Hence, if a receiver with coordinates (x,., z7) is located withincell ij, then define the dimensionless quantities p = (z? — x)/h and q = (z,. — z)/h. Theinterpolated traveltime is given byt(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, thenequation (8.8) is an exact expression for the traveltime at the interior point (x,., zr). In othercases, the interpolator (8.8) has accuracy 0(h2) (Dahlquist and Björk, 1974, p. 319) andthus is consistent with the level of accuracy associated with the forward modeling scheme.1758.4 Ray generationRaypaths are generated by following the steepest descent direction through a computedtraveltime field from each receiver back to the corresponding source point. This strategywas originally suggested by Vidale (1988), and has recently been implemented by Podvinand Lecomte (1991).The horizontal and vertical components of the traveltime gradient vector within cell ijare approximated by the centered finite-difference formulaeOt (t+1,+ t+,+i) — (tq + (8.9a)2h(t,1 + — (t +t1,) (8 9b2hAssignment of a constant traveltime gradient to a cell is compatible with the assumptionof locally plane wavefronts used in the forward modeling algorithm. The steepest descentdirection 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 raypathsthat cross cell ij have this same orientation angle. The lengths of the raypath segmentswithin 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 thecell at point A on its right boundary with coordinates (ma, za). Depending on the value ofthe steepest descent angle assigned to the cell, the ray may exit on any of the remainingthree sides or one of the four corners. The logic that selects one of these seven possibilities176xi xiii x—I Izizj+1zFig. 8.1. Raypath (heavy line) traced through square cell ij of a 2D slowness model. Rayenters cell at point A, follows the local steepest descent direction across the cell, and exitsat point B.is given in the first column of Table 8.1, where al and a2 are positive acute angles definedby—1 fZj+1 — Za\ —1 fZa Zj= tan h a2= tan hThese two angles are illustrated in Figure 8.1. After an exit option is selected, the coordinatesof 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 (columns4 and 5 in Table 8.1). In the particular case displayed, the raypath enters cell i,j + 1 onits top boundary. Hence, a different logical scheme is required to extend the raypath acrossthis next cell. A total of eight logical tables are necessary to handle all of the possibilities.B177Angle Range Exit Coordinates Next Cell Indiceso < ir/2 i + 1 j + 1ir/2 <8k <r — Xa + (Zj+i — za) cot 8k Zj+l i=lr—al Xi zi+i i—i j+11<9k.<lr+a2 xi zahtallOk i—i jxi Zj i—i f—ir + <8k <3ii-/2 x + (Zj — za) cot 0k i j — 13ir/2 9j <2ir z1 i + 1 j — 1Table 8.1. Ray tracing logic for an entry point on the right side of a square grid cell. Column1 gives possible ranges for the steepest descent angle 9k assigned to the cell. Angles al anda are defined in the text and ifiustrated in Fig. 8.1. Columns 2 and 3 give the horizontaland vertical coordinates of the ray exit point, respectively. (xa, za) are the ray entry pointcoordinates. 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, andthe 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 andexit coordinates are known. Although tomographic inversion requires only the value of thesegment length, the actual coordinates (xa, Za) and (zb, zb) are also retained in order tocreate raypath plots for diagnostic purposes.Ray tracing is initiated at each receiver position. If a receiver is located on a gridnode or grid line, then an average of the steepest descent angles from the neighboring cellsdetermines which cell the raypath enters first. Tracing then begins using the logic outlinedabove. However, if a receiver is positioned within a cell, a separate logical scheme generatesthe initial path segment to the enclosing cell boundary. The technique is similiar to thatdescribed above and is given in Appendix G. Iterative generation of path segments continuesuntil the raypath arrives at the boundary of a defined near-source zone. The size and shape178of this region vary slightly depending on whether the point source resides on a grid node, agrid line, or within a cell (see Appendix G). The final portion of the raypath is then taken tobe a straight line from the boundary point directly to the source position (c3,z3), regardlessof the local values of the steepest descent angle. This ray termination procedure is designedto overcome difficulties associated with nonuniformity of the traveltime gradient vector inclose 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 arrivalwavefronts 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 thehigh velocity zone. In addition, Figure 8.2b displays the raypaths traced through this timefield from 18 downhole receivers back to the source. The raypaths are orthogonal to thewavefronts, as expected.8.5 Inversion mathematicsAs 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 vector with relatively large and unrealistic cell-to-cell variations in slowness. The conditioningof these equations can be improved by incorporating model constraint information into theinversion. Hence, the system (8.3) is augmented with sets of linear equality constraints onthe updated model m + zm. In addition to mathematically stabilizing the inversion, theseequality constraints allow the convenient introduction of a priori geological or geophysicalknowledge (or bias) into the problem.In general, a set of linear constraint equations applied to the improved model m + mis written as the matrix/vector multiplicationB(m + z.m) = b, (8.11)179250horizontal position (m)Fig. 8.2. (a) Velocity model (contour interval = 100 m/s) with a horizontal low velocityzone 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 surfacesource. The raypaths are traced from 18 downhole receivers back to the surface source.000u00250 500rci)00tb 500180where the coefficient matrix B and the right hand side vector b are prescribed. Thesubscript m (n = 1,2,. .. , N) is used to refer to a particular set of constraints. A perturbationAm that simultaneously satisfies, in the least squares sense, the linearized data equationsand 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 therelative importance of the various terms. Extremizing with respect to Am yields thelinear algebraic equations[ATA + p, BBn] Am ATAt + — B m). (8.13)For nonzero ji, the coefficient matrix in this expression is usually nonsingular. The required model perturbation Am can be obtained by solving (8.13) using standard techniquesof numerical linear algebra. However, the coefficient matrix is large (K x K, where K is thenumber of slowness cells) and may be dense even if the original raypath and constraint matrices are sparse. Hence, it is advantageous to seek a solution method that avoids explicitlyforming the square matrices ATA andIt is straightforward to demonstrate that (8.13) are the normal equations associated withthe least squares solution of the rectangular systemA At—Bim)Am = . (8.14)vhi(bN- BNm)181Algorithm LSQR (Paige and Saunders, 1982) is used to solve equation (8.14) directly form. LSQR is an iterative solution technique for large and sparse systems of linear equationsthat is closely related to the conjugate gradient method. It is designed to seek the minimumnorm least squares solution of a set of equations. Numerical studies of LSQR applied tothe 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 matrixmay be large and sparse, a significant reduction in storage space is achieved by storing onlythe 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 particularlysimple to implement. -Finally, the correction to the slowness value at grid point ij is determined by averagingthe 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 modelare 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 nextiteration of the inversion can proceed.Several investigators apply a spatial filter to the slowness field between tomographiciterations in order to suppress short wavelength variations in the computed values (Dinesand Lytle, 1979; Radcliffe et al., 1984; Gersztenkorn and Scales, 1988). These variationsmay arise from noise in the traveltime data and/or instabilities in the inversion. Smoothingthen serves to condition the slowness field for the next forward modeling step. Hence, anoptional 9-point square smoother with prescribed weights Wi is included here. The smoother182is applied to the gridded values. The smoothed slowness value assigned to grid point ij isgiven bySj = Wi Si—1,j—1 + W2 Si,j_1 + W3 Si+1,j_IW4 Si—1,j + W5 S:j + W6 Si1,jtV7 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 conceptuaJly extending the grid by one point with the local values. In the examples to follow, theifiter weights are W5 4/12 and all other wj = 1/12.8.6 Model constraint equationsThere is a wide variety of linear equality constraints that can be employed in the tomographic inversion problem. The particular type of constraint used in this study is characterized 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 slownessmodel. Using the row-ordered indexing scheme, the kth component of the ifitered image Bmis given byCl 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 definedregion with the local cell slowness values. For example, application of the 5-point operatorto the upper left corner cell (i = 1, j = 1 corresponding to k = 1) of the model m is via theformulaci m2 + C2 ml + (c3 + C4 + c5)ml.One such constraint is applied to each cell of the slowness model. Development of the matrixrepresentation B for this set of equations is merely a matter of proper row and column183zizj+1 -zC3 C5 CiC2Fig. 8.3. Schematic representation of the 5-point constraint operator applied to cell ij of aslowness model. The slowness values in five neighboring cells are multiplied by constants Clthrough C5 and the products are summed.indexing. In general, the components of the right hand side vector b may differ, implyingspatial variation in the value of the applied model constraints.The 5-point operator is an extremely simple and flexible mechanism for introducingvarious 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 itsimmediate neighbors. Rather, inversion produces a model that is closest, in the leastsquares 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 finitedifference approximation to the horizontal derivative of the slowness field. A similarapproximation to the vertical derivative is achieved by setting c2 = 1/2h, C4 = —1/2hxi xi+i xI IC4184and all other ci = 0. These operators introduce first-difference regularization, or flattening, into the inversion.Case 3: ci = 1/h2, C5 = —2/h2,c = 1/h2 and all other cj = 0. This operator, as well as itsvertical counterpart, introduces second-difference regularization, or smoothing, into theinversion.Case 4: c = C2 = C3 = C4 = 1/h2 and C5 = —4/h2. This operator incorporates Laplaciansmoothing into the inversion.These constraints are inherently local in character; they consist of linear relations betweenadjacent 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 operatoroffers a reasonable compromise among the competing issues of flexibility, accuracy, and easeof implementation. -Combinations of constraints applied by the simple 5-point operator can also be considered. For example, the next section illustrates situations where both horizontal andvertical first-difference regulañzation is applied (i.e., N 2). The first system of constraintequations corresponds to c = 1/2h and C3 = —1/2h, while the second set corresponds toC2 = 1/2h and C4 = —1/2h. In each case, the right hand side vector b is set equal to 0.8.7 Synthetic examplesThe 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 mx 500 m velocity model used to generate synthetic traveltime data. Since the grid intervalis h = 5 m, there are IJ = 10201 grid points used for the forward modeling and K = 10000square cells used in the inversion.The data acquisition geometry used for the first example simulates a double-well VSPplus 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 m185separation) from each source position. In addition, the 9 downhole receivers in the right wellrecord traveltimes from 9 sources symmetrically placed in the left well. The complete setof 243 raypaths linking all source-receiver pairs is illustrated in Figure 8.4. Note that thesefirst arrival raypaths tend to avoid the low velocity zone, resulting in a region of reducedray coverage. Furthermore, no raypaths penetrate below 450 m depth.Contoured velocity tomograms obtained by inverting the combined VSP and crosswelltraveltimes are displayed in Figure 8.5. As indicated above, both horizontal and verticalfirst-difference constraints are imposed in the iterative inversions. Additionally, in Figure8.5b, the slowness values of the cells adjacent to the two boreholes are constrained by thetrue slowness function. In an actual field experiment, this information may be availablefrom borehole velocity logs. Each of these reconstructions produces an rms traveltirne errorof “-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 inversionsdepicted in Figure 8.5. The location and amplitude of the dipping high velocity anomaly areapproximately correct. Also, the shallow low velocity zone has been detected and correctlypositioned, although the actual velocity value at its center is about 90 rn/s too high. Thiseffect 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. Theprincipal difference between the two reconstructions appears in the region below 400 mdepth, where raypath coverage is negligible or nonexistent.Six (Figure 8.5a) and nine (Figure 8.5b) iterations are required to reduce the initialrms traveltime misfit to ‘—‘ 0.5 ms, with most of the improvement actually occurring on thefirst iteration. As stated above, the conjugate gradient solver LSQR is also an iterativealgorithm. Theoretically, LSQR requires at most K iterations to converge to the solutionof a system with K unknowns (assuming exact arithmetic can be performed). For thetomographic inversions described in this section, an acceptable solution is obtained withabout 10 iterations in LSQR, which is three orders of magnitude less than the theoreticalvalue K = 10000. This results in an appreciable saving in computation time.186Fig. 8.4. 243 raypaths traced through the velocity field of Fig. 8.2a. These rays link allsource-receiver pairs of a combined double-well VSP and crosswell experiment.0•4)O2a,C0LOO 250horizontal position (m)500187500250 500horizontal position (m)Fig. 8.5. Reconstructed velocity models (contour interval = 100 m/s) obtained by invertingthe 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 atboreholes. Max velocity 2801 m/s; mm velocity = 1724 rn/s.000tooa25000rhorizontal position (m)Fig. 8.6. Reconstructed velocity model (contour interval = 100 m/s) obtained by invertingnoise contaminated VSP and crosswell traveltimes. Max velocity = 2810 m/s; nun velocity= 1731 rn/s.188Figure 8.6 illustrates that the tornographic inversion is stable when the synthetic trayeltimes are contaminated with small amounts of random noise. Random numbers drawnfrom a uniform probability distribution on ±4 ms are added to the exact traveltimes. Theiterative inversion is initiated with the same constant slowness model used for the previous example. The convergence criterion for terminating iterations is arbitrarily selected tobe 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 theweights of the flattening constraints (..jij and in equation (8.14)) are set sufficientlyhigh.C00250 500189The final example examines the ability of the constrained inversion algorithm to reconstruct the interwell velocity field when only crosshole traveltimes are available. Figure 8.7displays the 81 crosswell raypaths. The zone of reduced ray coverage at shallow depths ismuch more extensive. Images obtained by inverting error free traveltimes are illustrated inFigure 8.8. Once again, both horizontal and vertical flattening constraints are applied. Theimage in Figure 8.8a displays a shallow low velocity zone and a deeper, dipping high velocityzone. However, the peak of the high velocity anomaly is shifted downdip by a significantdistance. An inversion including borehole velocity constraints is displayed in Figure 8.8b.- Surprisingly, the reconstruction is not improved; rather, a spurious region with high velocitygradients develops near the right well.The tomographic inversions in this section are all initiated with uniform slowness modelscalculated via the method described in Appendix H. Similar results are obtained from arange of nearby slowness values. Alternately, a nonuniform starting-mode1 obtained byhorizontal interpolation of the borehole slowness values can be be used. For the examplesconsidered here, this latter technique does not yield any obvious improvement in the velocityreconstructions.8.8 ConclusionFinite-difference traveltime computation offers an attractive alternative to conventionalraytracing for tomographic inversion purposes. The method is sufficiently rapid and accurate, and handles all of the various wave types that constitute first arrivals. Moreover, sincetraveltimes are computed throughout a slowness model, very general recording geometriesare easily accommodated. The main limitation of the technique is that it is restricted to firstarriving waves. Hence, the present formulation cannot be applied to the reflection tomography problem. However, current efforts to generalize finite-difference computation methodsto reflection traveltimes (e.g., Podvin and Lecomte, 1991) are encouraging.The introduction of constraint information into traveltime tomography is a responsibleway to address the nonuniqueness inherent in this inverse problem. Constraining informationmay arise from known geological or geophysical properties of the subsurface velocity modelDoqI-’. oodepth(rn)0 N 0 00tI) I-’.0C.. oC,’000 . 0 0 Ct,50025000191500Fig. 8.8. Reconstructed velocity models (contour interval = 100 m/s) obtained by invertingonly 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. Maxvelocity = 2862 m/s; mm velocity = 1728 rn/s.0-00050025000tO0 250horizontal position (m)192(say, from outcrops or borehole logs). Alternately, constraints may derive from a desireto impose certain reasonable attributes, like flatness or smoothness, on the constructedmodel. The method of incorporating constraints into the mathematical inversion procedureis adaptable to either viewpoint. Linear equality constraints are applied directly to theconstructed model, rather than to a model perturbation, and are satisfied in the least squaressense. The examples illustrate the imposition of flattening constraints and known boreholeinformation in the reconstruction of a smoothly varying interwell velocity field. Inclusion ofthis constraint information allows the solution of a problem that is strongly underdetermined- (243 data and 10000 unknowns). Nevertheless, a superior result is not necessarily achievedby the addition of the borehole velocity constraints.193CHAPTER 9SUMMARYThis 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 dimensional layered earth models in Chapters 2, 3, and 4. This class of models is characterizedby uniform velocity layers bounded by plane interfaces with arbitrary strike and dip. Computational procedures are developed for both single-layer and multilayered models. In thesingle-layer case, closed form mathematical solutions to the forward and inverse problemsexist. A rigorous derivation of a traveltime equation for critically refracted waves propagating 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. Thisformula forms the basis of an iterative head wave traveltime inversion algorithm designedto recover the parameters defining a single-layer or multilayered earth model. Inclusion ofconstraint information in the procedure, in the form of inequality relations satisfied by themodel parameters, often governs the ability of the algorithm to converge to a realistic solution. Tests with simulated and field data, acquired in various recording geometries, indicatethat the single-layer version of the algorithm is reasonably robust. However, the multilayered inversion algorithm appears to require fairly restrictive constraints in order to operateeffectively.(2) Improvements to various two dimensional refraction traveltime inversion methods thatallow 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 wavefrontmethod) 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 explicitlyincorporates undulating interfaces and horizontally varying refractor velocity into the model.Moreover, point values of interface depth, interface dip, and refractor velocity are obtainedfrom the observed arrival times. Hence, a depth profile of the critically refracting horizon can194be constructed via interpolation techniques. A rational procedure for calculating a smoothdepth profile, based on methods of linear inverse theory, is presented. This latter techniqueis not limited for use solely with the critical offset refraction profiling method, but mayaugment any traveltime analysis procedure that yields point depth or dip estimates of thesubsurface interface (such as the modified GRM developed here). Finally, an automatedimplementation of the classical wavefront method for interpreting refraction arrival timesindicates 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 velocitystructure. A simple imaging condition involving the reciprocal time then defines the locusof the subsurface refracting horizon. Tests with synthetic data indicate that both antidinaland synclinal structures can be imaged accurately. Shallow refraction data acquired at anarcheological site is also used to assess the workability of the algorithm.(3) Application of finite-difference traveltime computation methods to the two dimensionaltomographic inverse problem is treated in Chapter 8. An iterative algorithm is developed forreconstructing a. P-wave velocity field from measured first arrival times. Rapid and accurateforward modeling of all first arrival types (direct waves, head waves, and diffractions) forarbitrary source-receiver geometries is achieved with a finite-difference algorithm. Curvedraypaths, needed for converting traveltime residuals into localized updates to the slownessmodel, are generated by following the steepest descent direction through the computed trayeltime field from each receiver back to the source. Incorporation of constraint informationinto the procedure, in the form of horizontal and vertical first difference regularization, servesto 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 crosswellarrival times indicate that the algorithm can successfully reconstruct a smoothly varyinginterwell velocity field, even when the problem is severely underdetermined.These contributions can be applied to the solution of practical problems in seismicrefraction exploration, for various objectives and at many different scales. The stage is195now set for generalizing the techniques to more complicated three dimensional models withnonpiane interfaces and/or nonuniform velocities.196REFERENCESAckermanu, H.D., Pankratz, L.W., and Dansereau, D., 1986, Resolution of ambiguities ofseismic refraction traveltime curves: Geophysics, 51, 223-235.Adachi, R., 1954, On a proof of fundamental formula concerning refraction method of geophysical 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 forconstructing 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 Geophysicalprospecting: Am. Inst. Mm. Metallurg. Eng., 572-624.Berni, A.J., and Roever, W.L., 1989, Field array performance: theoretical study of spatiallycorrelated variations in amplitude coupling and static shift and case study in the ParisBasin: Geophysics, 54, 451-459.Berry, M.J,, and West, G.F., 1966, An interpretation of the first-arrival data of the LakeSuperior 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 velocityand depth in laterally varying media: Geophysics, 50, 903-923.Bois, P., La Porte, M., Lavergne, M., and Thomas, G., 1972, Well-to-well seismic measurements: 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 constantvelocity 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 ofstructure and velocity by seismic tomography: Geophysics, 51, 1559-1571.197Chiu, S.K.L., and Stewart, R.R., 1987, Tomographic determination of three-dimensionalseismic velocity structure using well logs, vertical seismic profiles, and surface seismicdata: Geophysics, 52, 1085-1098.Clayton, R.W., and McMechan, G.A., 1981, Inversion of refraction data by wave field continuation: Geophysics, 46, 860-868. VCunningham, 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 helpof 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 inversionof 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 seismicmethod, 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 emergedand 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-trimmedmeans: Geophys. J., 92, 67-72.Gjøystdal, H., and Ursin, B., 1981, Inversion of reflection times in three dimensions: Geophysics, 46, 972-983.Hadjidaki, E., 1988, Preliminary report of excavations at the harbor of Phalasarna in WestCrete: Am. J. Archaeol., 92, 463-479.Hagedoorn, J.G., 1959, The plus-minus method of interpreting seismic refraction sections:Geophys. Prosp., 7, 158-182.198Hales, 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 inversion: J. Can. Soc. Expi. Geophys., 20, 40-54.Hasselström, B., 1969, Water prospecting and rock-investigation by the seismic refractionmethod: 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: Geophysics, 47, 1431-1436.llatlierly, P.J., and Neville, M.J., 1986, Experience with the generalized reciprocal method ofseismic refraction interpretation for shallow engineering site investigations: Geophysics,51, 255-265.Hawkins, L.V., 1961, The reciprocal method of routine shallow seismic refraction investigations: Geophysics, 26, 806-819.Heiland, C.A., 1929, Modern instruments and methods of seismic prospecting, in Geophysical 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-dimensionalplane layer case: Geophysics, 41, 233-242.Hunter, J.A., and Pullan, S.E., 1990, A vertical array method for shallow seismic refractionsurveying of the sea floor: Geophysics, 55, 92-96.Johnson, S.H., 1976, Interpretation of split-spread refraction data in terms of plane dippinglayers: Geophysics, 41, 418-424.Johnson, L.E., and Gilbert, F., 1972, Inversion and inference for teleseismic ray data, inBolt, 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 refraction data: Bull., Seis. Soc. Am., 75, 865-880.Kilty, K.T., Norris, R.A., McLamore, W.R., Hennon, K.P., and Euge, K., 1986, Seismicrefraction 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 shallowtargets into the 1990s: Geophysics, 54, 1513-1542.199Lankston, R.W., 1990, High-resolution refraction seismic data. acquisition and interpretation,in Ward, S.H., Ed., Geotecimical and environmental geophysics, volume 1: review andtutorial: Soc. Expi. Geophys., 45-73.Lankston, R.W., and Lankston, M.M., 1986, Obtaining multilayer reciprocal times throughphantoniing: 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 undergroundimage reconstruction: IEEE Trans. Geosci. Remote Sensing, 18, 234-240.Macrides, C.G., Kanasewich, E.R., and Bliaratha, S., 1988, Multiborehole seismic imagingin steam injection heavy oil recovery projects: Geophysics, 53, 65-75.McMechan, G.A., Harris, J.M., and Anderson, L.M., 1987, Cross-hole tomography forstrongly 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 problemof 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 refractionmethod: 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., Seismictomography with applications in global seismology and exploration geophysics: D. ReidelPubi. Co., 1-23.Ocola, L.C., 1972, A nonlinear least-squares method for seismic refraction mapping - partII: 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 firstarrivals: Geophys. Prosp., 37, 455-465.Paige, C.C., and Saunders, M.A., 1982, LSQR: an algorithm for sparse linear equations andsparse least squares: ACM Trans. Math. Software, 8, 43-71.Palmer, D., 1980, The generalized reciprocal method of seismic refraction interpretation:Soc. Expi. Geophys.200Palmer, D., 1981, An introduction to the generalized reciprocal method of seismic refractioninterpretation: 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 shallowrefraction seismology: Expi. Geophys., 21, 33-44.Palmer, D., 1991, The resolution of narrow low-velocity zones with the generalized reciprocalmethod: Geophys. Prosp., 39, 1031-1060.Parker, R.L, Shure, L., and Hildebrand, J.A., 1987, The application of inverse theory toseamount magnetism: Rev. Geophys., 25, 17-40.Pavlenkin, A.D., Saiculina, S., and Vinnik, A.A., 1986, Method of conjugate points in theinterpretation of traveltime curves of refracted waves: Izvestia, Earth Physics, 22, 795-803.Phadke, S., and Kanasewich, E.R., 1990, Seismic tomography to obtain velocity gradientsand three-dimensional structure and its application to reflection data on VancouverIsland: Can. J. Earth Sci., 27, 104-116.Podvin, P., and Lecomte, I., 1991, Finite difference computation of traveltimes in very contrasted 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 techniquefor 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 wavesin 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 refraction 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 inverseproblems: 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 theinterpretation of data from seismic surveys: Geophysics, 22, 9-21.201Schenck, 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 oldmine workings by the seismic refraction method, in Ward., S.H., Ed., Geotechnical andenvironmental 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 incomplex 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 lineprofile 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: Geophysics, 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 ofmultilayer 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 Archregion: Ph.D. thesis, Univ. of British Columbia.Zelt, B.C., Effis, R.M., and Clowes, R.M., 1990, SCoRE ‘89; the Southern Cordillera refraction experiment; description of data set: Lithoprobe report, 20, Univ. of BritishColumbia.202APPENDICESAppendix A: Conic sections in poiar coordinatesExpressions for conic sections in plane polar coordinates can be found in most textbookson analytic geometry. However, these formulae typically assume that the coordinate originis 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. Theprinciple axes of the conic are rotated by an arbitrary angle S with respect to the coordinateaxes, 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 theunprimed system through the angle 9. The coordinates of a given point relative to eachsystem 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 sexuimajoraxis 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; inthe case, the effipse encompasses the origin. In the primed reference frame, the equation ofthe 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.203However 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 xX cos a and y = X sin a. After some algebraic reduction, a quadratic form in the radialcoordinate 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 azimuthalangle 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 asspecial cases of equation (A3). Thus, setting l = 0 (i.e., origin coincident with the ellipsecenter) immediately yieldsI 1—e2X(a) = a41V l—ecos(a—0)Also, if l = a e, then the origin is located at a focus of the ellipse. In this case (A3) gives1—e2X(a)=This form is common in mechanics; it describes the trajectory of a particle moving underthe influence of an inverse square central force.A similiar analysis reveals that equation (A3) also applies to a hyperbola with the samecenter location and principle axis orientation. The only difference is that the numerical valueof the eccentricity exceeds one (e > 1) for a hyperbola.204xx,0yy,Fig. Al. A noncentered, rotated effipse.205Appendix B: Traveltime analysis for a simple 3D modelA compact formula for head wave traveltime in a simple three dimensional earth modelis derived in Chapter 2 (equation (2.10)). The purpose of this appendix is to verify that themore general expression (3.12) reduces to this known result in the special case of a two mediamodel 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 inFigure 2.2) is confined to a single plane. The two vectors n2 and p22 form an orthonormalbasis for all vectors in this plane. Hence, the propagation vectors p12 and q12 can be resolvedalong this basis. Geometric analysis of Figure 2.2 yieldsP12 (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 byrecognizing thatrQ — rpP22 =trQ—rp1I206where rp and r are the position vectors of points P and Q in Figure 2.2. Thusrq — 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)n2P22 = . (B4)cosSubstituting these expressions for 2 and P22 into equations (B2) yields the ratherformidable formulaesin #2 COS 82 cos(i + 8) + sin i cosP12 1cos S+ sin #2 sin 82 cos(i + 8) + sinicsin’11 + [cos #2 cos(i + 8) k, (B5)cosS j L cosSand—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 cosSIt 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 that207P12,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 resultT2(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 traveltimeformula (2.10) previously developed for this simple model situation.208Appendix C: Traveltime analysis for a simple 2D modelA useful check on the validity of the general head wave traveltime formulae in Chapter 3 isprovided by analyzing a specific situation for which a closed form traveltime solution exists.This is especially important because several of these expressions disagree with analogousformulae published recently by Diebold (1987). The inconsistency between the two resultsis evident when source and receiver are located on separate interfaces with different dipangles. In this case, it is not possible to effect a coordinate frame rotation such that bothinterfaces become horizontal (i.e., parallel to the zy plane).Consider the simple two dimensional earth model depicted in Figure Cl. Two subsurfaceinterfaces have the same dip angle ç. The source S is located at the coordinate origin on thesurface and the receiver R is located on subsurface interface 2. L is the source-receiver rangemeasured parallel to the dipping interfaces. The perpendicular distances from the origin tothe subsurface interfaces are d2 and d3; vertical layer thicknesses at the same point are h1and h2.The traveltime of a head wave formed on interface 3 can be derived from first principlesof refraction traveltime analysis. Rotating the perspective through the small angle cp transforms the situation into a one dimensional problem. The total traveltime is easily obtainedby summing the refraction delay time associated with each layer that the wave traverses,together with L/v3. The result isT3 =. +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 ço209/Fig. Cl. A simple 2D earth model with three media. S is a surface source and R is a buriedreceiver.where ZR is the horizontal coordinate of the receiver. Substituting these expressions intoequation (Cl) yieldsT3(ZR) = ZR +h1 cos(013— +2h cos 823 COSV3COS9 1)1 V2(C2)The traveltime predicted by the general formula (3.16) is now compared with this specifictraveltime solution. For the situation being examined, (3.16) simplifies considerably. Thefollowing conditions hold:Sd1XR/d2viI//RV3(p210(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 yieldsT3(R) = XRCOS 50 q23, + SIR 50 q23,z +h1 P13,z+h2(p3,—q23,z) (C3)V2COSSO V2Since the earth model is strictly two dimensional, and the recording proffle is oriented normalto the common strike direction of the interfaces, the unit propagation direction vectors arecontained entirely within the xz plane. Thus, the components of these vectors can beobtained 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) givesZR Sm 623 h1 cos(613— ) 2h cos COS 50T3(ZR)= + +v2cosço v1 V2But since sin 923 = v2/v3, this equation reduces immediately to= ZR+h1 cos(913—0)+2h COS 023 COS (C4)vacOsço vi V2211This is in complete agreement with the known solution (C2)! In contrast, equation (21) inDiebold (1987), when applied to the model depicted in Figure Cl, becomesT3(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 correcttraveltime only if the dip angle ço equals zero, i.e., all interfaces are horizontal.212Appendix D: Derivation of formula (6.13)A standard technique from the calculus of variations is used to demonstrate the vadidityof equation (6.13). If the function m”(z) is perturbed by an arbitrary, but small, amountSm”(x), then a variation S is induced in the functional 4(m”). An expression for thisvariation is obtained by evaluating equation (6.12) at m” + Sm”. Hence5cJ (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 smallperturbation, this can be achieved only ifm”(x) =---[cobs — erd(m”)] WTWP(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 differencebetween observed and predicted data. Parker et al. (1987) obtain an analogous result for adifferent inverse problem formulation.213Appendix E: Vanishing of two derivativesThis appendix proves that the extremum of the objective function in equation (6.17) isinvariant with respect to the two abscissae c and C2. Differentiating 4’ with respect to theseparameters yields the two equations= a’ [2r WTW + I] aad 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 ijc9d2+ 2 [cilh(b) — d2 [k(b) — (b — ci)h(b)] — eobs]Ta, (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 vectora: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)214T’C2 C2Expressions for the partial derivatives of the inner product matrix F are required. Differentiating 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 derivativesof the kernel function vector p(x; ci, c2) are obtained by differentiating expression (6.10).Hencec2)= 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, substituting (E8) into the analogue of (E6) for ÔF/0c2 gives the result0C2= 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)215Relations (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 parametersd1 and d2, then aTh(b) = aTk(b) = 0. Thus, equations (Eli) and (E12) reduce to thedesired results0, — 0. (E13)Note that both of the conditions a Th(b) = 0 and a Tk(b) = 0 are necessary for 8c1/8 2 tovanish.216Appendix F: Spatially correlated traveltime errorsSuppose that x is a vector of ii random variables with covariance matrix C. Then, thevector of n random variables y given by the linear transformation y = Ax possesses thecovariance matrixC1, = ACXAT. (Fl)If C,, and C1, are prescribed, then equation (Fl) can be solved for the (n x n) transformationmatrix A.Since the covariance matrices are symmetric and positive definite, they may be factoredvia Cholesky decomposition as follows:i _y r—.LJ 12Z ‘-‘1, — U1, .LJ1,,where L,, and L1, are lower triangular matrices. Then, one solution of (Fl) isA = L1,;, (F2)which can be readily verified by substitution.An important special case occurs when the z1 are a set of independent random variablesdrawn from the one dimensional normal distribution with zero mean and unit standarddeviation. In this situation, the covariance matrix C,, equals the identity matrix. Thetransformed random variables are given simply by y = L1, x. Moreover, each y is normaflydistributed 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 ar.ndom time error with an assigned standard deviation o. A double-tailed exponentialfunction is used to evaluate the correlation coefficient between the arrival time errors at217geophone stations i and j (Berni and Roever, 1989). The elements of the covariance matrixbecome{cgj = exp(—d/D)ojoj, (F3)where djj is the distance between stations i and j, and D is an adjustable parameter calledthe correlation distance. A large value for D (relative to the dj) implies that the individualarrival 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 withzero mean and unit variance. Hence, the y calculated via the above procedure become a setof correlated, normally distributed traveltime errors.2i8Appendix G: Ray initiation and termination logicIf 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 descentangles 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 containedin Table Ci is used to determine the initial path segment and the cell that the raypathsubsequently enters. Angles /9i through /38 are all positive acute angles and are illustratedin Figure Gib.Now assume that points 1, 2, and 3 are sources, rather than receivers. Then, theset of immediately adjacent cells constitutes a near source zone. Tracing of an incomingray continues until the boundary of this zone is encountered. For example, the four cellssurrounding point 1 in Figure Cia form a square with side length 2h. For source points 2and 3, the near source zones are rectangles with vertical and horizonal dimensions of 2h x hand h x 2h, respectively. Finally, if the source is located within a cell (Figure Gib), theniterative ray tracing proceeds to the boundary of the enclosing square cell with side lengthh. In all cases, the final raypath segment is obtained by drawing a straight line from theboundary point directly to the source position.ZihZrZj+13 ah1 .2hh h219Fig. Gi (a) Points 1, 2, and 3 denote receivers (or sources) located on a grid node or gridlines. (b) A receiver located within a square grid cell.X1 Xr xi+1h220Angle Range Exit Coordinates Next Cell Indiceso <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+1Bk=lr/+/3 Zj+1 i—i j+11—/4<Bk<+/5 Z zT+(zj—xf)tanOk i—i jSk=1+/35 Z Zj i—i j—131r/2—/6<8k <31r/2+/37 xr+(zj—z,.)cotOk zi I j—1Sk=31r/2+/ i+1 j—i2ir—/38 <Bk <2w xj+1 z,. + (zj± — zr)tan 0k i + 1 jTable Gi. Ray initiation logic for a receiver located within a cell. Column 1 gives possibleranges for the steepest descent direction 6k assigned to the cell. Angles /9j through /38 areillustrated in Fig. Gib. Columns 2 and 3 give the horizontal and vertical coordinates ofthe 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.221Appendix H: Optimum starting slownessIn the absence of a priori information, a. uniform slowness model is used to initiate thetomographic inversion procedure. The value of this slowness can be chosen to minimizethe rms traveltime residual on the first iteration. This tends to reduce the total number ofiterations required for convergence to a specified misfit level.Let t be the th observed traveltime, and let the straight line distance between theassociated source and receiver be dj. Then, the difference between the observed traveltimeand the time predicted by a uniform slowness s is t — s dj. The 12 norm (squared) of allsuch differences is4(s) = (td)2where N is the total number of observed arrival times. Extremizing 4(s) with respect to 8yields-- EfL1jts—st—rAlternately, the problem can be posed in terms of an optimum velocity instead of anoptimum slowness. The result is v = 1/s. For the combined VSP and crosswell datasetsused 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 Issued | 1992 |
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 |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-12-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052957 |
URI | http://hdl.handle.net/2429/3268 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1992-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1992_spring_aldridge_david_franklin.pdf [ 3.87MB ]
- Metadata
- JSON: 831-1.0052957.json
- JSON-LD: 831-1.0052957-ld.json
- RDF/XML (Pretty): 831-1.0052957-rdf.xml
- RDF/JSON: 831-1.0052957-rdf.json
- Turtle: 831-1.0052957-turtle.txt
- N-Triples: 831-1.0052957-rdf-ntriples.txt
- Original Record: 831-1.0052957-source.json
- Full Text
- 831-1.0052957-fulltext.txt
- Citation
- 831-1.0052957.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052957/manifest