GEOMETRICAL FORMULATION AND TRAVELTIME COMPUTATIONIN HETEROGENEOUS LAYERED MEDIAbyXINGONG LIB.Sc., Tongji University, P. R. China, 1984A THESIS SUBMI TIED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE FACULTY OF GRADUATE STUDIESDepartment of Geophysics and AstronomyWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIADecember 1992©Xingong Li, 1992In 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 of6- E-0 p H y SICS^A S Ton n/o/v)/ The University of British ColumbiaVancouver, CanadaDate TAN01R1 2 , 0'93DE-6 (2/88)ABSTRACTThis thesis deals with the subject of seismic traveltime computation. The velocity field isrepresented by an arbitrary distribution of rectangular or square cells. The cells with an interfacecrossed are represented by two non-rectangular cells. Each cell is assigned a constant velocity.A geometrical approach is presented to construct the Vidale formula (Vidale, 1988) and the lineartraveltime interpolation, LTI, formulae (Matsuoka et al., 1990; Asakawa and Kawanaka, 1991;Podvin and Lecomte, 1991; Schneider et al., 1992). The geometrical construction allows an easyvisualization of the principles which underlie both the Vidale and the LTI methods. From thisgeometrical point of view, subdividing the cell boundary into smaller segments (the subdividedLTI, or the SLTI method) improves the approximation of the plane wave assumption and hence,improves the resulting accuracy.A new approach is developed to calculate the traveltimes for a heterogeneous model withinterfaces. The algorithm involves three steps:1. The discretization of the velocity field using square or rectangular cells. For cells on aninterface, each cell is divided into two secondary non-rectangular cells. The two non-rectangular cells are assigned different velocities, so that even an interface with a largecontrast in velocities can be modelled.2. The calculation of traveltimes at the interface locations.3. The computation of the reflection traveltimes and the transmission traveltimes by expandingthe wave fronts from the interface points upwards and downwards from the interface.The principles used in this approach are Fermat's principle and Huygens' principle. Theknown grid points or interface points are treated as secondary sources and the traveltimes atthe points within a cell are calculated and revised by finding the shortest traveltimes fromthe known points. Because of the special consideration given to the interface cells, the newapproach is found to be more accurate than the Vidale and the LTI approaches. With respect toiireflection traveltime computation, the new algorithm calculates the reflection traveltimes for allthe receiving points at the same time and is, consequently, faster than the reciprocity principlemethod (Matsuoka and Ezaka, 1991; Podvin and Lecomte 1991).Various examples using the algorithm are presented. Model experiments show that thecomputation errors are smaller in a high velocity field than in a low velocity field. The rectangulargrid is tested by increasing the vertical sampling rate of the square grid by a factor of two. Theresults show that computational errors decrease in the horizontal direction and increase in thevertical direction.CONTENTS^ iABSTRACT^ iiLIST OF FIGURES ^ viACKNOWLEDGEMENTS viiiChapter 1^INTRODUCTION^ 1Chapter 2 METHODS OF TRAVELTIME COMPUTATION ^ 42.1^The Huygens' Principle Method^ 42.2 The Vidale Method ^ 42.3^The Linear Traveltime Interpolation (LTI) Method ^ 5Chapter 3 GEOMETRICAL FORMULATION^ 93.1^The Vidale Formula^ 93.2 The Non-subdivided LTI (NLTI) Formula ^ 93.3^Relationship between the Vidale and NLTI Formulae ^ 123.4^The Subdivided LTI (SLIT) Formula ^ 12Chapter 4 TRAVELTIME COMPUTATION 15Chapter 5 EXAMPLES ^ 245.1^Model 1 — Uniform Half Space ^ 245.2^Model 2 — Vertical Gradient Velocity Model^ 245.3^Model 3 — Dipping Interface ^ 335.4^Reflection Traveltime Computation 365.5^Model 4 — Traveltime Computation for Converted Waves ^ 405.6^Effects of the Magnitude of the Velocity Field ^ 44iv5.7^Rectangular cells ^ 455.8 Model 5 — Heterogeneous Layered Medium ^ 475.9^Discussion ^ 50Chapter 6^CONCLUSIONS ^ 516.1^Geometrical formulation 516.2^Traveltime Computation in Heterogeneous Layered Media ^ 516.3^Reflection time computation ^ 526.4^Effects on the Computation 536.5^Application in Traveltime Tomography^ 53REFERENCES ^ 54Appendix A Possible Interface-cell Combinations ^ 57Appendix B Traveltime Computation for a Point in a Vertically Increasing VelocityModel ^ 58LIST OF FIGURES1. The Vidale method ^ 62. The LTI method 63. The geometrical formulation of Vidale's formula^ 104. The geometrical formulation of NLTI formula 115. The geometrical formulation of SLTI formula ^ 146. Different types of irregular interface cells 167. Two examples of interface cell decomposition ^ 168. Wave front expansion in the presence of an interface 229. Triangular cell traveltime minimization ^ 2310. Traveltime contours for model 1 2511. Traveltime errors for model 1 at x = 0 ^ 2712. Traveltime contours for model 2, 50m x 50m grid ^ 2813. Traveltime errors for model 2 at x = 0, 50m x 50m grid 2914. Traveltime contours for model 2, 5m x 5m grid ^ 3115. Traveltime errors for model 2 at x = 0, 5m x 5m grid 3216. Traveltime contours for model 3 ^ 3417. VSP traveltimes for model 3 at x = 0 3518. Reflection traveltime contours for model 3 ^ 3819. Reflection traveltimes and their errors for model 3 at the surface ^ 3920. Reflection and converted traveltime contours for model 4^ 4121. Reflected and converted traveltimes and their errors for model 4 at the surface ^ 4222. VSP traveltime errors for model 4 at x = 0^ 4323. VSP traveltime errors for model 1 at x = 0 with different velocities ^ 4424. Traveltime contours for model 1 with rectangular grids ^ 46vi25. Traveltime computations for rectangular cells^ 4726. The velocity field of model 5 ^ 4827. Traveltime contours for model 5 4828. Traveltimes at surface z = 0 and at borehole x = 0 for model 5^ 49viiACKNOWLEDGEMENTSI am very thankful to my supervisor Dr. Tadeusz J. Ulrych for his guidance, encouragementand the support to this study. Special thanks are due to him for the time and the patience devotedin helping the writing. I wish to thank Toshifumi Matsuoka, Eiichi Asakawa, David Aldridgeand John Hole for generous and stimulating discussions concerning all aspects of traveltimecomputation. A traveltime computation program using Vidale's method was provided by DavidAldridge and the LTI traveltime computation program used for model 1 was provided by EiichiAsakawa. I thank Yao-guo Li for his help with the diagrams. I am thankful to the wholedepartment of Geophysics and Astronomy for all the help I have received during my stay. Also,thanks to all my friends in Canada. Your friendship makes the country even more beautiful.Thanks to my parents, my grandfather, my sisters and all my friends in China for all thoseletters with so much love and friendship.viiiChapter 1INTRODUCTIONThis study can be divided into two parts. In the first part, a geometrical approach is presentedto formulate the two currently used finite difference, FD, methods — the Vidale and the LTImethods. In the second part, a new algorithm is developed to calculate the traveltimes in anarbitrarily distributed velocity field with the presence of interfaces.Traditionally, traveltime computation for an arbitrary velocity distribution requires someform of raytracing during the computation. Two basic raytracing approaches are the shootingmethod and the bending method (e.g.,terven$, et aL, 1977; Thurber and Ellsworth, 1980). Inthe shooting method, an initial value problem is solved. A ray is shot from a starting point withan initial take-off angle, then a fan of rays is generated by changing the take-off angle. Thecorrect ray path and traveltime may be approached when the end point is close to the receivingpoint. The bending method starts the raytracing with an initially estimated ray path betweena source and a receiver. Subsequently, the ray path is bent by a perturbation method until itsatisfies a minimum traveltime criterion.Some difficulties exist in these raytracing methods:1. For most of surface seismic records, vertical seismic profiles (VSPs), or cross-hole seismicrecords, the traveltimes that can be easily picked are those of the first arrival events, therefore,first arrival time computation is of important in traveltime tomography. However, forstrongly varying velocity fields, many ray paths between a source and a receiver usingraytracing methods may be generated and hence the first arrival times may be missed.2. When there are many source points and each source point has many receivers, the raytracingmethods can be computationally inefficient.123. For very complicated velocity structures, it may be difficult to find the ray path. A smallchange in the take-off angle may result in a large change in the ray path, and it may takemany trial and error attempts to find the correct ray path.To overcome some of the difficulties mentioned above, three kinds of traveltime computationapproaches have been developed. The first approach is the Huygens' principle method (Itoh andYamada, 1983) which treats each grid point as a secondary source. At a certain grid point, allpossible traveltimes between this grid point and the secondary sources around it within a certainrange are considered. The shortest traveltime is selected as the correct traveltime. Saito andOhtomo (1986) and Saito (1989,1990) implemented this method using graph theory and Moser(1991) attacked the problem by means of a network technique. Vidale (1988) presented anothertraveltime computation method which solved the eikonal equation using a FD approach. VanTrier and Syme (1991) also used a FD approach in their upwinding method. Matsuoka et al.(1990), Asakawa and Kawanaka (1991), Podvin and Lecomte (1991) and Schneider et al. (1992)developed a method of computing the shortest traveltime by determining the minimum ray pathover an individual cell. Matsuoka et al. (1990) and Asakawa and Kawanaka (1991) namedthis approach the linear traveltime interpolation, LTI, method. These authors also extendedthe LTI approach by subdividing each cell boundary into several segments which improves theaccuracy of the traveltime calculation. This method is called SLTI in this thesis. The formulaused by Podvin and Lecomte (1991) and Schneider et al. (1992) is a special case of the oneused by Matsuoka et al. (1990), Asakawa and Kawanaka (1991), without the subdivision ofthe cell boundary. This non-subdivided LTI is called NLTI in this thesis. The Vidale andthe LTI methods are based on the plane wave assumption. Special techniques are required forspherical wave front expansion which will not be considered in the research. In these methods,the velocity field, or alternatively, the slowness field, is represented by many discrete rectangularor triangular cells with a constant velocity or slowness in each cell. One of the difficulties in3the cell method is the problem of representing the velocity field when interfaces exist. When acell is crossed by an interface with a high velocity contrast between the two sides, the cell cannot be properly represented by a constant cell velocity. In this thesis, non-rectangular cells areused to represent the interface cells. A new algorithm is developed to calculate the traveltimesin a layered earth model.A geometrical formulation, which is a new method to derive the Vidale formula and the LTIformulae, is introduced in this thesis. A clear understanding may be obtained of both Vidale'sand the LTI formulation by this means. The geometrical point of view also allows a deeperinsight into the traveltime computation and helps in simplifying and optimizing the traveltimecomputation algorithms.The reflection traveltime computation using the reciprocity principle method has been studiedby Matsuoka and Ezaka (1991) and Podvin and Lecomte (1991). Two difficulties exist inthis method. First, for a shot point, the reflection traveltimes for the receiving points arecalculated separately. This is inefficient when many receivers are considered. Second, thereis no proper handling of the interface cells. The reflection traveltime accuracy may suffer fromthe misrepresentation of the interface cell velocity and the interface position. The new algorithmhas overcome these problems.This thesis is arranged as follows. In chapter 2, the three traveltime computation methods,i.e., the Huygens' principle method, the Vidale method and the LTI methods, are briefly reviewed.In chapter 3, the Vidale formula and the NLTI and SLTI formulae are derived using thegeometrical formulation method. A numerical example is used to show that subdiving thecell boundary in the LTI method makes the computation closer to the plane wave assumption,thus improving the computing accuracy. Chapter 4 introduces the new algorithm for traveltimecomputation. Finally, in chapter 5, five models are used to test the new algorithm.Chapter 2METHODS OF TRAVELTIME COMPUTATION2.1 The Huygens' Principle MethodAn approximate approach using Huygens' principle regards every grid point as a secondarysource and the remaining grid points as receiving points. It takes into consideration all of thepossible ray paths between the source grid point and the receiving point and, at each grid point,chooses the shortest time from all the source points as the correct traveltime. This method iseasy to implement and can produce a good approximation to the ray paths in a relatively complexslowness field. The simplest eight direction pattern is tested in chapter 5. The Huygens' principlemethod is a fast approach, but the errors are too large for an eight direction basic pattern. Theaccuracy of the method can be improved by increasing the number of directions of the basicpattern (Saito, 1989). The problem is that the basic pattern with more than eight directionscan not deal with cell velocity changes within the covered range of the basic pattern. Thebasic patterns are also affected by the grid types. Mandal (1992) recently applied the Huygens'principle method to triangular gridded models.2.2 The Vidale MethodVidale (1988) solved the eikonal equation•j; =S2(X7Z)(12 (aty(1)4using a square grid with an average slowness. As illustrated in Figure 1, ti , t2, t3, t4 aretraveltimes at the four corners, ti , t2, t3 are known and t4 is the required traveltime. s isthe average slowness and h is the length of each side. The derivatives are expressed asat^t4+t2- ti -t,ax 2hat^t4 + t3 - ti - t2•az 2hSubstituting equation (2) into equation (1), we obtain Vidale's formula14 = (sh)2 — (t3 — t2)2 -Vidale expanded the wave fronts from small to larger traveltimes as illustrated in Figure 1(b). Heimplemented the algorithm in three steps. First, the relative minima are found for the calculatedgrid points (filled circles). In this simple case, only one minimum exists. The traveltimes of thegrid points just outboard of those at the relative minima (the filled circles at the outer ring inFigure 1(b)) are calculated using a noncentered finite-difference of equation (1). The plane-waveformula he used is t3 = to + V(sh)2 — 0.25 (12 — 11)2, where to is the relative minimum timein the inside ring, t3 is the time to be found, and t1 and t2 are the times on either side ofthe point whose time is to For the next two steps, the traveltimes at the hollow circles arecalculated using equation (3) in both directions leaving the minimum points. For example, forthe bottom edge of the outer ring, Vidale solved the unknown times rightwards until either arelative maximum or the edge is encountered. Then, the unknown times are solved leftwardsuntil reaching a relative maximum or edge. The difficulty is that when the velocity contrastbecomes large, equation (3) is violated and causality problems arise. This aspect of the Vidalemethod has been studied by Qin et al. (1990, 1992).23 The Linear Traveltime Interpolation (LTI) MethodPodvin and Lecomte (1991) showed that there is another way of using the eikonal equationto compute the traveltime t4 at the rectangular grid point C, as shown in Figure 2(a), from thecomputed traveltimes ti and 12 at grid points A and B. Suppose that the mesh slowness is s5(2)(3)Etox-- IA t B t 2raytb6• time calculated(a)^0 time not calculated^(b)Figure 1 Vidale's method uses (a) traveltimes of three grid points to calculatethe traveltime of the fourth grid point, and (b) square wave front expansion.Dt31'4-- a(a) C t4(b)Figure 2 (a) LTI approach minimizes t4 as a function of x to find the traveltimeof the third grid point. (b) Except the ray path considered in (a), other ray pathshave to be considered. There are 16 possible ray paths to a grid point. Nine of theseare plotted here. The minimum of the possible traveltimes is the correct traveltime.7and that the ray goes through the cell boundary AB to point C. The partial derivatives can beexpressed as = 1.4 t2 . Here x is the horizontal axis and z is the vertical axis.Combining the above expressions with the eikonal equation, we havet4 = t2^Vs2a2 (t2 ii ) 2 .a(4)Equation (4) is also derived by Schneider et al. (1992). Their dynamic programmingtraveltime computation which is based upon Fermat's principle, minimizes a function of theray position over one side of the cell. Suppose, initially, that the ray goes through side AB atpoint E with time to which is assumed to be the linear interpolation of ti and t2. Referringto Figure 2(a), we obtaint2 —to = tl ^x.aThe traveltime at grid point C can then be written as a function of the position x ast4 = to + .51/b2 — (a — x)2t2 — t1t1^X + Vb2 — (a — x) 2 .aMinimizing (6) as a function of x we haveb(t2Vs2a2 — (t2 — ti) 2Substituting equation (7) into equation (6), we obtain equation (4). The correct point on thecell boundary AB should be the one at which the ray path satisfies Snell's law when it crossesthe boundary. From this point the ray reaches point C and its traveltime, according to Fermat'sprinciple, is the minimum. For plane waves, the linear interpolation assumption, i.e., equation(5), becomes true and, therefore, the correct time is obtained. Since the wave front can beapproximately but not exactly a plane wave during the computation, the computed time canx = a(5)(6)(7)8only be an approximation of the correct time, and the computed time is always greater thanthe true time.Clearly, the ray may also cross side AD or travel along the cell boundary DC or BC.There are, in fact, as pointed out by Podvin and Lecomte (1991) and illustrated in Figure2(b), 16 possibilities. These authors have developed a massively parallel and a sequentialimplementation of this approach.Matsuoka et al. (1990) and Asakawa and Kawanaka (1991) have divided each cell boundaryinto several smaller segments and have assumed that the traveltime at any point on the cellboundary can be linearly interpolated between traveltimes at adjacent segment points on the samecell boundary. Their method is called Linear Traveltime Interpolation (LTI) method because ofthis assumption. Obviously, the approach of Schneider et al. can be regarded as a special case ofLTI without the subdivision of the cell boundary. We will call the non-subdivided LTI formulaNLTL and the subdivided LTI formula SLTI. The SLTI formula is given later in equation (20).Non linear interpolation has been suggested by Schneider et al. (1992).Matsuoka et al. (1990) and Asakawa and Kawanaka (1991) have demonstrated that with thesame number of cells the SLTI approach has a higher accuracy than the Vidale method. Thismeans that SLTI can achieve the required accuracy using larger cells, therefore fewer modelparameters. As is well known, fewer model parameters can save computer space and computingtime in the inverse problem. When interfaces exist, some difficulties occur in representing themodel using square or rectangular cells. We will introduce irregular cells, which can further bedecomposed into triangles and rectangles, to represent the discrete model.Chapter 3GEOMETRICAL FORMULATIONThis chapter presents a geometrical formulation of equations (3) and (4), assuming planewave fronts. The subdivided LIT approach (SLTI) can also be derived without any difficulty.As will be shown, this geometrical view of the ray tracing formulation can give us a betterunderstanding of the plane wave assumption. Geometrical approach has been used by Dix(1981) in 1952 to formulate the eikonal equation.3.1 The Vidale FormulaSuppose that ti, 12, t3, t4 are the traveltimes at the four corners of a cell. The plane wavefronts associated with a square cell are shown in Figure 3. With h the length of each side ands the cell slowness, we have that EC = (14 — ti)ls is the interval between the wave front at14 and ti, and BF = (13 — 12)Is is the interval between the wave front at 13 and 12. SinceAC I BD and ABDF = AACE, we obtain(14 -4 2 (t3 - t2)22= (V2h) . (8)Equation (8) leads immediately to Vidale's formula expressed by equation (3).3.2 The Non-subdivided LTI (NLTI) FormulaFigure 4(a) shows the geometry associated with a rectangular cell. a and b are the lengthsof the two sides. Initially we assume that the ray traverses side AB to reach the grid pointC. The LTI formula uses 1 and 12 to calculate 13 and can be easily deduced using the910Figure 3 The geometrical formulation of Vidale's formula.geometrical approach. In Figure 4(a), we see that BE = (12 — ti)/s, FB = (14 — 12)/sand AABE N LBCF, where "—" stands for similarity.Consequently, CF/BE = b/a and CF = b(12 — t i )/as. From ABCF, we obtain(14 _ 12) 2 b2 ( b ) 2 02 ti )2(9)s )^a)^s )Rearranging equation (9) leads to the LTI formula given by equation (4).So far, we have supposed that the ray goes through cell boundary AB. From the geometricalpoint of the view we can easily deduce which side of the cell boundary the ray actually traverses.For example, a critical case is showed in Figure 4(b). Referring to this figure, we see that, inabAABE, AE = /a2 —^)2. Further, AE • AC = ab and consequently, AEs 12 -1-b2 •We now have the following possibilities:1. 12 — t =^, when the ray passes through A;2. t 2 — ti^a1::■■142 , when the ray passes to the right of point A;3. t2 — ti > :772 sEb2 , when the ray passes to the left of point A;4. The ray does not passes through the cell in which case we must use neighboring cellsto calculate 14.1 1(a)^ (b)Figure 4 (a) Geometrical formulation of the NLTI formula. (b)A critical situation — ray passes through point A to point C.Knowing the traveltimes for any two points in the cell, the geometrical approach, assumingplane wave fronts, allows us to compute the traveltimes for any other points inside the cell oron the cell boundaries. For example, referring to Figure 4(a), using t1 and t3, we can calculatet4 as followsHC = (t3 — ti) ,bsC4 —st3) 2a2 (HC)2We can combine the above two equations to obtainat4 = t3 —b Vs2012 (t 3 — 4)2'We can also compute t3 from ti and t2 by means ofb It3 = +^S2a2 — (t2 — )2.(10)(12)1233 Relationship between the Vidale and NLTI FormulaeAs a further example, suppose 12 and 13 are known and we wish to compute 14. SinceABCF ACDH, we haveCFs b_ .t4— t3^aFrom ABCF,(CF)2 (14 — t2 ) 2 = b2.Combining the above two equations yields(1 + e2) ti — 2 (e2t3 t2) t4 e2t3^— b2s2 = o^ (15)where e = bla. For a square cell, a = b = h,e = 1, and12 t3V2h2s2= 2^— (12 — t3)2 .2 (16)Assumingt4^= t2 t3)^ (17)equation (16) leads to Vidale's formula, i.e., equation (3).3.4 The Subdivided LTI (SLTI) FormulaMatsuoka et al. (1990) and Asakawa and Kawanaka (1991) subdivided each cell boundaryinto several segments in order to improve the accuracy of the NLTI method. Without thesubdivision, i.e., the formulae which we deduced in the previous section, the NLTI method isless accurate than the Vidale approach. With several subdivisions however, the SLT7 version ofthe LTI method achieves a higher accuracy than Vidale's method. This is shown with a simpleexample in Figure 5(a). Assume the cell size is of unit length and width, and that the slowness(13)(14)d V^s2d2 _ (tA — iB )2Iz -F z tp — tA (18)13is 1.0. The units are ignored here for convenience. Assume that the traveltimes at the gridpoints and the segment points along the cell boundaries EF and FC are known and considerthe calculation of the traveltime at grid point D. The Vidale method uses the traveltimes at gridpoints E, F and C. NLTI uses the traveltimes at grid points E and F and the three segmentSLTI uses the traveltimes at segment points A and B. The calculated traveltimes at grid pointD are 2.287, 2.324 and 2.246 respectively. The errors relative to the true traveltime, 2.236, are,0.051, 0.088 and 0.010 respectively. From the geometry of the ray path, OD, we see that whenwe subdivide the cell boundary into three segments the plane wave front assumption is bettersatisfied within each individual segment than before subdivision, so that a higher accuracy isachieved. We conclude that the subdivision of the cell boundaries better satisfies the plane waveassumption. From Figure 5(b) for the plane wave case, we can deduce the SLTI formula usingthe geometrical reconstruction without any difficulty. Since L\ABK AGDH,Also /ABK N RAGE, andto — tBlx— = \s2 _72a^(tA — tB) 2Combining the above two equations we obtain Asakawa and Kawanaka's formula4^lx (4^4 \ 1 Z^.,12^(4^4 N2tp to + (tA — GB) +^2 — (tA — GB) .(19)(20)14d^lxI -4—o- 1-4-1.- 1 lz(a)^ (b)Figure 5 The calculation of traveltime at grid point D. (a) Circular wave front case.Vidale's formula uses the traveltimes at grid points E,F and C. LTI uses thetraveltimes at grid points E and F without subdivision or segment points A and Bwith three subdivided segments for each cell boundary. (b) The plane wave frontcase. SLTI formulae can be deduced from a geometrical consideration (see text).Chapter 4TRAVELTIME COMPUTATIONAn excellent example of the application of FD traveltime calculations to tomographicinversion is the recent work of Aldridge and Oldenburg (1992). These authors used smoothestmodel constraints to solve the underdetermined inversion problem. Although such a constraint iscertainly reasonable for the cells within layers, it is not applicable to the layer interfaces wheresharp velocity contrasts may exist. As seismic sections indicate, in most cases of interest, theearth may be adequately represented as a layered structure. In order, therefore, to represent theslowness field in the real earth more closely, one can either decrease the cell size or use anirregular mesh. Decreasing the cell size will obviously increase the computation time of bothforward and inverse calculations and could result in a considerably increased computational load.Reflection traveltimes have been calculated using the reciprocity method of Matsuokaand Ezaka (1990, 1992) and Podvin and Lecomte (1991). Unfortunately, because of theimproper representation of the velocity field along the interface, these methods confront thesame difficulties on the interface as do the transmitted traveltime computations.To overcome the above problems, a new algorithm is developed with the help of a geometricalunderstanding of the traveltime computation formulae. The new algorithm allows the expansionof the wave front using irregular cells and thus can overcome the above problems without theneed of finer sampling. The algorithm consists of three steps. First, the model is discretizedinto rectangular or square cells. A cell crossed by an interface is subdivided into sub-cells withdifferent slownesses. The second step is started in the layer where the source is located. Thetraveltimes of both the grid points and the interface points in this layer are calculated in thisstep. A modified NLTI sequential method is applied for traveltime computations in rectangular1544444.411%%116.4463/4444,4interfaceFigure 6 Different types of irregular cells when an interface exists. Allthe possible interface and cell combinations are considered in AppendixA. Details corresponding to the shadowed cells are shown in Figure 7.interface16(a) (b)Figure 7 Two typical irregular cells. We can decompose these cells into (a)one rectangular and one triangular cell, both with slowness s 1 and anothertriangular and rectangular pair with slowness s 2 , (b) three rectangular andone triangular cell with slowness s i and one triangular cell with slowness s 2 .cells and non-rectangular cells. In the third step, the interface grid points with correspondingtimes already calculated in the second step, are treated as secondary sources. The reflection,refraction, converted wave and the transmission times are calculated by expanding the wavefronts from the interface points. The steps of the algorithm will now be described in detail, anda few examples are presented in the next chapter.17Stepl. The interface discretizationThis algorithm assumes a layered earth model. The model is divided into / columnsand J rows, with K = I J cells. The kth cell can also be identified by its upperleft grid point (i, j). The relationship between the cell number and the grid index isk = (j —1)- I-1-i. The slownesses inside the layers are represented by the discretizedcell slownesses. The interface between two layers is represented by many straight linesegments as shown in Figure 6. When an interface crosses a cell, the cell is dividedinto two smaller non-rectangular cells with two different slownesses. These kinds ofirregular cells can be further decomposed into secondary rectangular and triangularcells as shown in Figure 7. All the possible interface segment positions in a cell areshown in Appendix A. An array is used to describe the cell types. Rectangular cellsand non-rectangular cells are treated differently and the approach is described in thefollowing steps.Step2. From source to the interfacesa) Rectangular cellsInitially, the shot point traveltime is assigned time zero and the rest of the gridpoints are assigned a large number. Podvin and Lecomte's (1991) massively parallelor sequential implementation is readily applicable to guarantee the least traveltimesfor the grid points in the source layer. The sequential implementation is chosen tominimize computational expense. The sequential implementation expands the wavefronts from the shot point to all the grid points inside the current layer and on theinterface. The expansion is stopped when the interfaces or the model boundaries arereached. The traveltime computations in a rectangular cell are described here and thecomputations in the interface cells are detailed in the next paragraph. Figure 8 choosesa particular row of cells to show how the wave fronts are expanded. It is supposedthat the expansion is in the downward direction and the source is somewhere at the18top of this row of cells. The traveltimes of the upper grid points , filled circles, havealready been calculated, and the traveltimes of the bottom grid points and the interfacegrids, hollow circles, are to be calculated. In order to satisfy Fermat's principle,the traveltimes of the grid points in the cells have to be minimized after the wavefront expansion. The algorithm consists of both leftward and rightward computationswhich are done cell by cell, as shown in Figure 8. The reason for using leftward andrightward computations is to expand the wave fronts from the minimum time pointsand thereby satisfying the causality condition. In the leftward computation, the wavefronts are expanded from the minimum time points to the left, and in the rightwardcomputation, the wave fronts are expanded from the minimum time points to the right.Each computation requires two steps to guarantee that each grid point has the leasttime. For example in the upper part of Figure 8, during the leftward minimization incell MDAB where M is the grid point with minimum time, the two-step approachis as follows:1. The minimization of traveltimes at point B and D. The direct times from M to B andD are computed and compared with the previous times. If the former are shorter thanthe latter then the times are updated by the times from M to guarantee B and D havethe least times. The ray paths considered are labeled by the number "1" in Figure 8.This step can be expressed by the equationt = min(t, t min(s) • h),where t is the time of either B or D, tm is the time at point M, min(s) gives thesmaller slowness between the current cell and the adjacent cell of either side MD orside MB, and h is the corresponding length of the cell.2. The minimization of the time at grid point A. The possible ray paths are the direct raypaths from point D, M and B, and the ray paths from side DM and MB using the19NLTI formula. The traveltimes of the above possible ray paths are computed and theleast time is saved.b) Irregular cellsOn the left side of Figure 8, the interface cuts the cells into two non-rectangularcells with different slownesses. The traveltimes of the irregular cells in the first layerare to be minimized. To do so, an irregular cell is first decomposed into rectangularand/or triangular cells as shown in Figure 7(b). Suppose that ti, t2 and t3 are calculatedand t2 is the least time. The upper pentagonal cell is decomposed into three rectangularcells and one triangular cell. Next, t5 and to are calculated by linear interpolation ofthe adjacent grid points. Then, the traveltimes of four secondary cells are minimized.The sequence of the secondary cell traveltime minimization depends on which pointhas the minimum time. In this example, if t2 is the minimum time, the minimizationshould be in the sequence of the upper-right rectangle first, the other two rectanglessecond, and the triangle last. If t3 or t7 is the minimum, the sequence should bethe lower-right rectangle first, the triangle second, and the other two rectangles last.If ti or t8 is the minimum, the sequence should be the upper-left rectangle first, thetriangle second, and the other two rectangles last. The rectangular cell minimizationis described above.Figure 9 shows how to minimize the traveltimes of a triangular cell. First theminimum time is found. If t7 or t8 is the minimum time, a two-step approach similarto the one adapted for a rectangular cell has to be applied. The ray paths consideredin the two steps are labeled "1" and "2". For example, in Figure 9(a), if t8 is theminimum, the two steps are as follows.1. At the to point, the direct traveltime from t8 point is computed and compared with theprevious to. The smaller time is kept. The same computation is repeated for the t7 point202. The possible traveltimes to the 17 point are calculated and compared with the previous17. The least time is kept. The possible traveltime taken from the diagonal side to theto point is also computed and is compared with the previous time and the smaller timeis kept.Similar considerations are required if 1 7 is the minimum time, as shown in Figure9(b). If to is the minimum as shown in Figure 9(c), only the direct ray paths need beconsidered.Step3. From the interface to the neighboring layersThe interface points are treated as secondary sources in this step and the transmis-sion, reflection and converted wave traveltime are then computed by expanding thewave fronts from the interface points.1. To calculate the transmitted traveltimes, the wave front expansion is carried on from theinterface points to the grid points in the next layer until the next interface is reached.2. To calculate the reflected traveltimes back to the present layer, a large number is setfor all grid points with the exception of the interface points and the wave fronts areexpanded from the interface upwards.3. To calculate the converted wave traveltimes, the slownesses of the corresponding wavehave to be substituted into the computation.The wave front expansion has to be started from the minimum time points on theinterface in order to obviate the causality problem. Hence, the minimum times at theinterface are found first and the expansions are started from the minimum points.Step4. Multi-layer traveltime computationThe above method can be easily applied to multi-layer media. First, the traveltimesof the points at the closest-to-source interfaces are calculated using the approachdeveloped above. Then the traveltimes at the remaining grid points can be calculated21by expanding the wave fronts from those interface points to the remaining interfacepoints. The further that the interface is from the source the better the approximationto the plane wave assumption and hence the smaller the computational errors.• Traveltime computed0 Traveltime to be computedFigure 8 Wavefront expansion when an interface exists.22 23t 8It 7(a) t 8 minimumt 7(b) t 7 minimumt 7(C) to minimumFigure 9 Triangular cell traveltime minimization.24Chapter 5EXAMPLES5.1 Model 1 — Uniform Half SpaceIn the following example an uniform half space model is used to compare the LTI, Vidale'sand the Huygens' principle method. In the case of an uniform slowness field, the cell slownessis the exact representation of the field and no discretization errors exist.The model used here is a 500m x 500m square region divided into 50m x 50m smallersquares. The source is located at (500, 50). The velocity is 2000m/s. The wave fronts of thecalculated traveltimes using the above three methods are plotted in Figure 10 with the dottedcurves representing the true wave fronts. In order to compare the accuracies of the methods, theerrors between the computed times, using different methods, and the true times are calculated atthe receiving points along the vertical line x = 0 and are plotted in Figure 11. Both the NLTIand SLTI versions of the LTI method are used here. The SLTI method uses three subdivisionsand has the best accuracy among the four results. As discussed before, the SLTI accuracy ishigher because SLTI satisfies the plane wave front assumption better than other methods. Themaximum errors with respect to the true times are computed as a measure of the accuracy. Thelargest errors are 4.83, 0.67,1.76, 22.16ms for NLTI, SLTI, Vidale's, and the Huygens' principlemethod, respectively. These correspond to 19.32%, 2.68%, 7.04%, 88.64% of 25ms which is thetime taken by a ray to traverse a cell in a direction parallel to a cell boundary.5.2 Model 2 — Vertical Gradient Velocity Model5.2.1 A 50m x50m gridThe true traveltime for a given grid point in an arbitrary distributed velocity field is difficult500100 150 200 250 300 350 400 450horizontal position (m)0501001502005 250300a.a)10 3504004505000NLTI and true time (dx=dz=50m)50 100 150 200 250 300 350 400 450 500horizontal position (m)(a)SLTI and true time (dx=dz=50m)25(b)Figure 10 Comparison of the FD traveltime contours for model1: (a) NLTI, (b) SLTI. The dotted lines are true traveltimes.Vidale's and true time (dx=dz=50m)o50100150200E 2506 300a.o.)1:1 35040045050 100 150 200 250 300 350 400 450 500horizontal position (m)(C)Huygens and true time (dx=dz=50m)50 100 150 200 250 300 350 400 450 SOOhorizontal position (m)(d)Figure 10 (cont.) (c) Vidale's method, (d) Huygens' principle method.26271,1.1,17,1.11.^SLTI11111■111/ \-W..--/^\^Via-tire's. //^. itt-4 \It^NLTI^.... . -. .Ale.....".., Huygens' principle,e`.......^-.-•a^I^1^■^■^I^.^I^1^a^1^1^1^1^I^I^1^i^,^I^I^■^1^I^I^I"0 100 200 300 400 500depth (m)Figure 11 Traveltime errors for model 1 at x=0.to formulate analytically, therefore the accuracies of the FD traveltime computation methods aredifficult to estimate. In this section we choose a constant vertical gradient velocity model tocompare the NLTI and the Vidale traveltimes with true traveltimes. The reason for using NLTIinstead of SLTI is for the sake of simplicity. The accuracy of Huygens' principle method canbe improved by considering more complex basic patterns with more directions (Saito, 1989),but these patterns can not handle the velocity changes along the newly added directions. Forthese reasons, only the NLTI and the Vidale methods are compared in later sections. The exacttraveltime formulae of any grid point are derived in Appendix B. The velocity of a grid pointis expressed as v(x, z) vo kz with vo the velocity on the surface, k the vertical velocitygradient. The discretized cell velocity is the average of the velocities at four grid points. Welet k = 4.0/s, and vo = 1800m/s, so for a 50m x 50m grid, the grid velocities are 1800m/sfor the grid points at z = Om, 2000m/s for the grid points at z = 50m and 3800rrzis for thegrid points at z = 500m. The cell velocities, therefore, are 1900m/s for the cells at the top,2100m/s for the cells at the second row and 3700m/s for the cells at the bottom row of the grid.The wave fronts calculated by NLTI and Vidale's method are plotted in Figure 12. Thetraveltime errors at x = 0 are plotted in Figure 13. In Figure 12(a), NLTI traveltimes showrelative large errors for the grid points at the top of the model around z = 50m. This phenomenon50 100 150 200 250 300 350 400 450 500horizontal position (m)NLTI28^I0^50 100 150 200 250 300 350 400 450 500horizontal position (m)(a)Vidale's method(b)Figure 12 Comparison of the FD traveltime contours for model 2: (a) NLTI, (b) Vidale'smethod. The dotted lines are true traveltimes, The grid size is same as Figure 10.290^100^200^300^400^500depth (m)Figure 13 Traveltime errors for model 2 at x=0, 50m x 50m grid.can be explained as follows. The velocity function is represented by the discretized cell velocitiesduring the NLTI computation. Since the velocity is assumed constant in each cell, it is only anapproximation to the continuous model. As mentioned above, the cell velocities are 1900m/s atthe first row, 2100m/s at the second row, so that the direct wave traveltimes will be calculatedusing the faster velocity, 2100m/s. In the near source region, these direct waves travel fasterthan the refracted waves and thus are the first arrivals. They also travel faster than the truedirect waves and thus the errors are positive.5.2.2 A 5m x 5m gridThe natural way to mitigate the discretization errors is to decrease the grid size. A denser5m x 5m grid has been tested and the results are plotted in Figures 14 and 15. The absolutevalues of the maximum errors are listed in the table below, where the percentages are withrespect to the cell transit time, i.e., 25ms for a 50m x 50m grid and 2.5ms for a 5m x 5m grid.It is obvious that decreasing the cell size can improve the traveltime accuracies for bothmethods. It seems that for this model, the improvement in the Vidale method is greater thanthat for the NLTI method. For both grids, the accuracy of the Vidale method is superior to thatof the NL11 method. This suggests that for the continuous model, the NLTI method can not30be improve on the Vidale approach. The same result holds for the constant velocity model asshown in the previous section.NLTI Vidale'sMax. error Percentage Max. error Percentage50 X 50 5.61 22.5% 2.76 11.0%5 X 5 0.68 27.2% 0.20 8.0%^.Table 1 Comparison of the maximum relative errors between NLTI andthe Vidale method. The first column indicates the cell size. Max. erroris the maximum relative error between the computed time and the truetime. Percentage is the Max. error with respect to the cell transit time.I^I^I^I^1^1___.50 100 150 200 250 300 350 400 450 500horizontal position (m)NLTI (dx=dz=5m)31I^1^I^4^1^I^I50 100 150 200 250 300 350 400 450 500horizontal position (m)(a)Vidale's method (dx=dz=5m)(b)Figure 14 Comparison of the traveltime contours formodel 2, same as Figure 12 but with grid size 5m x 5mo 100 200^300depth (m)400 500:.5v-0.6a)32Figure 15 Traveltime errors for model 2 at x:1 with grid size 5m x 5m3353 Model 3 — Dipping InterfaceUsing the algorithm developed in the previous chapter, the traveltimes are calculated for amodel with two uniform layers separated by a dipping interface. The geometry is the same asbefore. The interface is defined by two end points with coordinates (0, 77.5) and (500, 377.5),about a 30° dip. The cells are 5m x 5m. Since the new algorithm is a combination of the NLTIformula and the special treatment for interface cells, NLTI is still used to indicate the algorithm.Figure 16 shows the wave fronts computed by the NLTI and the Vidale approach and Figure17 shows the traveltime curves and the errors with respect to true times at the receiving pointsalong the vertical line x = 0. The traveltime from the source transmitting through an interfaceto a given point can be formulated analytically (Aldridge, 1992).There are three kinds of waves in Figure 16: the direct wave form the source, the headwave in left part of the first layer and the transmission wave in the second layer. As for thecomputation of direct wave traveltimes, the Vidale approach is more accurate than the NLTIapproach. For the head wave and transmission wave, because of the difficulties in handlingthe interface cells during the computations, the Vidale approach is less accurate than the NLTIapproach . Figure 17 shows that the largest errors along x = 0 are —0.77ms for NLTI and+1.77ms for the Vidale method which correspond to 30.8% and 70.8% of the cell transit timerespectively. The errors for the constant velocity model in the previous section are 27.2% and8.0%. A comparison of the above errors shows the effects of the interface. It is obvious thatthe existence of the interface does not affect the NLTI method too badly, but affects the Vidaleapproach significantly.NLTI (dx=dz=5m)(a)(b)Figure 16 Traveltime contours for model 3. The dotted lines represent the truetraveltime wave fronts. (a) NLTI and true times, (b) Vidale's and true times.34100 200^300depth (m)ig' 1oo0o$.4 oo......litv -1E:ao 100_fj true timesA Vidale method• LTI method^--.^i^,^.^.^1^I^I^I^I^i^I^I(a)200^300depth (m)(b)400 500500400Figure 17 Traveltimes at x=0 for model 3-- a vertical seismicprofiling example. (a) the traveltimes, (b) errors to the true times.35365.4 Reflection Traveltime ComputationWe have demonstrated the traveltime computation using the LTI and the Vidale methodsfor an uniform medium with constant velocity, a continuous medium with a vertical velocitygradient and a layered medium with a dipping interface. All the examples illustrate first arrivaltime computations. The first arrivals may consist of the direct, the transmitted, refracted ordiffracted waves. Because of the importance of the reflection data in exploration seismology, anexample of reflection traveltime computation is considered in this section. The forward reflectiontime computation is important in seismic reflection traveltime tomograghy and is also useful inwave- equation migration.Mtsuoka and Ezaka (1990, 1992), Podvin and Lecomte (1991) have developed the reciprocityprinciple method to calculate reflection traveltimes. The reciprocity principle method requiresthree steps to calculate the reflection traveltimes:1. the computation of the traveltimes from a source point to all grid points;2. the computation of the traveltimes from a receiving point which is set as a temporary sourcepoint to all grid points;3. the addition of the two traveltimes calculated in step 1 and step 2 at the interface points. Theminimum added time along the interface points is the reflection traveltime from the currentsource to the current receiver according to Fermat's principle.Two difficulties exist in this method. First, for the case that many receivers exist for eachsource, the reflection time for each receiver is calculated by repeating the last two steps andis, consequently, computationally expensive when many receivers are considered. Second, ifrectangular or square grids are used and the interface is defined by the cell boundaries, theinterface position may not be well defined, which leads to errors. A denser grid is requiredin order to achieve a higher accuracy but implies a greater number of model parameters andconsequently a greater computational load.37The traveltime computation approach developed in the previous chapters represents theinterface by grid points of irregular cells and expands wave fronts from source point to interfacegrids and then from the interface grids to all the receiving points. Consequently, this approachovercomes the shortcomings of the reciprocity principle method. Figure 18 shows the wavefronts calculated by the new method for model 3 which has been used in the previous section.Figure 19 shows the traveltimes at the surface and the corresponding errors with respect to truetraveltimes.The computed traveltimes represent reflection times on the right side of the model and head(refracted) wave times on the left side of the model. The refraction times are computed alongthe interface where the velocity of the second layer is used for the computation. If only thereflection times are required, the refraction time computations can be eliminated.The larger error at the left ending point is caused by the causality problem. Causality implies,as stated by Vidale (1988), that the time for the part of the ray path leading to a point mustbe known before the time of the point can be found. For this model, the refracted ray pathsactually require some grid points outside of the left model boundary. Since these grid pointsare not available, large errors occur.Except for the receiving point at the origin, the maximum reflection time error is 0.71mswhich is at the right end of the surface. This error is 28.4% of the transit time and is about thesame as the error in transmission traveltime computation. Figure 18 shows that the reflectionwave fronts are close to plane waves in the top layer, hence the error in the reflection traveltimecomputation is small. The head wave fronts are plane waves, which will not contribute errorsto the head wave traveltime computation. For the above reasons, accuracy on the surface isrelatively good.450 500100 150 200 250 300 350 400horizontal position (m)NLTI (dx=dz=5m)Figure 18 Reflection traveltime contours for model 3 byNLTI method. The same geometry as in Figure 16 is used.381.111111f^I300^ true timeNLTI[h) 2804h)260240^111.1.1121■11.111111.111.11■^(1100^200^300^400^503horizontal position (m)39(a).IIIIIII■1.11. 100^200^300^400^500horizontal position (m)(b)Figure 19 Reflection traveltimes of a common shot point gather for model 3at z=0. (a) The traveltimes and (b) the errors with respect to the true times.405.5 Model 4 — Traveltime Computation for Converted WavesP-SV converted waves are generated by the conversion of compressional energy into shearenergy at the interface and are often observed, for example, in three component vertical seismicprofiling Inversion of the shear velocity, Vs, is of importance in many different aspects ofseismic exploration. In order to perform the inversion, forward modelling is required. Thissection shows an example of the P-SV reflection traveltime computation for model 4.The model parameters in model 4 are the same as model 3 except that the velocities of thetwo layers are altered to vi = 6000m/s and v2 = 2000m/s. Since the high velocity layer is ontop of the lower velocity layer, there is no refraction on the interface in model 4. The Vp/Vsratio of the top layer is set to 2.0, so the shear velocity is 3000m/s.Figure 20 shows the P-SV and P wave traveltime contours and a comparison with the truetimes. Figure 21 shows the traveltime curves and the errors computed on the surface. Figure21(b) illustrates that the error of the P time and the error of the P-SV time are about the same.Some interesting differences occur when the velocities are altered. Comparing the reflectionerror of model 4 in Figure 21(b) with the reflection error of model 3 in the right part of Figure19(b), a significant improvement in the accuracy is seen. Comparing the transmission time errorof model 4 in Figure 22 with model 3 in Figure 17(b) shows a similar improvement. Furthertests to determine the effects of velocity are described in the next section.oso1001502005 250.4....1 30004a)ij 3504004505000^50 100 150 200 250 300 350 400 450 500horizontal position (m)(a)50 100 150 200 250 300 350 400 450 500horizontal position (m)(b)Figure 20 Traveltime contours for model 4. The dotted lines represent the true traveltimewave fronts. Traveltime contours in the first layer are: (a) P-SV wave and (b) P wave.41100^200^300^400^500horizontal position (m)(a)100^200^300^400^500horizontal position (m)(b)Figure 21 Surface traveltimes for model 4. (a) P and P-SVtraveltime curves, (b) errors with respect to the true times.424300 -0.11,4 -02"C)a) -0.3.0-OA0 500400100 200^300depth (m)Figure 22 VSP traveltime errors with respect to true traveltimes at x=0 for model4. Comparing with Figure 17(b) of model 3, the errors in this diagram are smaller.400 5000^100^200^30depth (m)o00 4CO14IDfa-4 -64i-o..-IioCD -80..41-10445.6 Effects of the Magnitude of the Velocity FieldBoth reflection and transmission time errors show the differences between model 4 andmodel 3. This indicates that the magnitude of the velocity field effects the accuracy of traveltimecomputation. In order to examine the relationship between the traveltime error and the velocity,several constant velocity models are used in this section. The errors with respect to the truetimes are computed and plotted in Figure 23. As can be seen, the error decreases non-linearlywhen the velocity increases.Figure 23 VSP time errors of model 1 with different velocities. From bottom totop the velocities change from 1000 m/s to 6000m/s with a 1000m/s increment.455.7 Rectangular cellsIn this section, 50m x 25m and 25m x 50m rectangular grids are used for traveltimecomputations in the constant velocity model. The results are compared with those of a 50m x 50msquare grid.Figure 24 shows the traveltime contours calculated using rectangular grids. The velocity is2000m/s which is the same as that of model 1. Comparisons with Figure 10(a) which illustratesthe results of square grid indicate some changes in the error distributions. For instance, theaccuracies are improved in upper left part of Figure 24(a) and are degraded in the remainingpart of the figure. These changes can be explained with the help of Figure 25.A 50m x 50m square cell is divided into two 50m x 25m rectangular cells in Figure 25.The unknown grid points are the hollowed circles and triangles. The times at the filled circlesand triangle are used by the NLTI formula to calculate the unknown times. In Figure 25(a), thetraveltimes at the shorter cell boundaries are used for the computations. The traveltimes are closeto plane waves at these shorter boundaries, so higher accuracies are achieved. This is similarto the SLTI method as discussed in section 3.4. For the grid point at the left-bottom cornerof Figure 25(b), the computation is carried out in two steps for the rectangular grid. First, thetraveltimes at two small triangles are calculated. Then the traveltimes of the two triangles areused to calculate the time of the left-bottom point. The two step computation is done using onestep in the case of a square grid. The error at the left-bottom corner is the sum of errors in thesetwo steps. The wave front approximation in each of the two steps is about the same as usingone square cell, so each step yields about the same error as using one square cell. Therefore theerror using a rectangular grid is about twice the error obtained using a square grid.(a)50 100 150 200 250 300 350 400 450 500horizontal position (m)50 100 150 200 250 300 350 400 450 500horizontal position (m)(b)Figure 24 Traveltime contours using rectangular grids. (a) 50m x 25m, (b)25m x 50m. A constant velocity model is used. The velocity is the same as in model 1.4647(a)^ (b)Figure 25 The accuracies of traveltime computations in rectangular cells.The dotted line cuts a square cell into two rectangular cells and adds twogrid points (triangles). (a) Using the times of the grid points at the short sideof the cell boundary achieves a higher accuracy. (b) Using the times of thegrid points at the longer side of the cell boundary achieves a lower accuracy.5.8 Model 5 — Heterogeneous Layered MediumThe final model considered is model 5 which is model 3 with three added velocity anomalies.The anomalies are defined by the coordinates and the amplitude values of the peaks. Theanomalies decay exponentially. The first two anomalies are in the first layer, with a 2000m/sbackground velocity, centered at (200,100) and (400, 200) with —1000m/s and 1000m/smaximum anomaly values respectively. The third anomaly is in the second layer, with 6000m/sbackground velocity, centered at (100,100) with —2000m/s in magnitude. The velocities in thefirst layer are changing from 1000m/s to 3000m/s, and the velocities in the second layer arechanging from 4000m/s to 6000m/s. Figure 26 shows the velocity field of the model. Figure27 illustrates the traveltime contours for model 5. The contours in the first layer correspond tofirst arrivals of reflection or refraction waves and the contours in the second layer are those ofthe transmitted wave. The traveltimes at the surface and at the left side borehole (x = 0) areplotted in Figure 28 as solid lines and are compared with the traveltimes of model 3 which arerepresented by dotted lines.t` 400020006 0a'Cr,Zi" ° ^ 0.....^"Cr,^o^0^0o 0 odepthFigure 26 Velocity field of model 5. The model is model 3 with three anomalies added.Figure 27 Traveltime contours for model 5.480 100^200^300^400^500horizontal position (m)400 5000^100^200^300depth (m)(a)(b)Figure 28 Traveltimes at (a) surface, z=0; (b) borehole, x=0. The solid lines representthe traveltimes of model 5, and the dotted lines represent the traveltimes of model 3.300250200240g 220(1.) 200918016049505.9 DiscussionWhen the traveltimes have been computed, the ray paths can be found either by following thesteepest gradient in traveltime from the receiver back to the source (Vidale, 1988; Aldridge andOldenburg, 1991) or by the backward computation method developed by Matsuoka et al. (1990)and Asakawa and Kawanaka (1991). Once the traveltimes and the ray paths are known, velocityperturbations about the model can be inverted using various techniques. The model is thenrecomputed using these perturbations and the final result is achieved after a number of iterations.The problem for tomography is that the interface structures have to be provided in advancein order to compute the traveltimes in the forward problem. One way to handle this problem isto combine reflection tomography with depth migration (Bording et al., 1987). Another way is toset the interface depths as model parameters and linearize the relationship between the interfacedepth perturbations and the traveltime residuals. This method can be used to invert the interfacedepths and the model velocities simultaneously (e.g., Bishop et al., 1985; White, 1989).Chapter 6CONCLUSIONS6.1 Geometrical formulationThree kinds of first arrival time computation methods in an arbitrarily distributed velocitymodel, the Huygens' principle method, the Vidale method, and the LTI (both NLTI and SLTI)method, have been briefly reviewed. All of these methods are based on the plane waveassumption. The Vidale formula can be derived by solving the eikonal equation using thediscrete FD method. The LTI formulae can be derived either by the discrete FD method orby the least time functional minimization method which guarantees that Fermat's principle issatisfied within each individual cell. A geometrical formulation approach has been presentedto derive both the Vidale formula and the LTI formulae. From the geometrical point of view,the reason that the SLT1 formula is more accurate than the NLTI and the Vidale formulae canbe clearly interpreted. For the SLTT method, dividing the cell boundary into several segmentsmakes the wave front at each segment closer to a plane wave, thus increasing the accuracy ofthe computation. A simple numerical example was used in section 3.4 to illustrate this fact.6.2 l'raveltime Computation in Heterogeneous Layered MediaWith the help of the geometrical understanding of the FD traveltime computation formulae,a new algorithm was developed to compute the traveltimes for a heterogeneous layered velocitymodel. The new algorithm includes three steps.511. Representation of the cells that cross an interface by two non-rectangular cells. The non-rectangular cells can then be decomposed into secondary rectangular cells and/or triangular52cells. A cell traveltime minimization method was developed to guarantee that the traveltimeat each secondary grid point is minimum, thus satisfying Fermat's principle.2. Computation of the traveltimes at the interface points.3. Calculation of the reflection times and the transmission times by expanding the wave frontupwards and downwards from the interface points. Converted traveltimes can be calculatedby substituting compressional wave velocities with shear wave velocities.The algorithm was tested for models with only one interface, but extension to multi-layeredmodels is straightforward. Multiple reflection times can then be calculated if the computationbetween two interfaces is carried out upwards and downwards several times. No subdivided cellboundaries were used in the algorithm for reason of simplicity.Instead of using non-rectangular cells only at the interface cells, step 1 can be extended touse non-rectangular (polygonal) cells for the entire model. This is suitable for a velocity fieldwhich can be represented by many blocks, each having a constant velocity. Each side of thepolygon should be short enough to adequately conform to the plane wave assumption. If anyside of the polygon is too long, it should be subdivided. The traveltime minimization methodfor triangular and rectangular cells can be similarly applied to polygonal cells.63 Reflection time computationA comparison with the reciprocity principle method showed that the new algorithm is fasterand more accurate. The algorithm can compute reflection times from a shot point to all thereceiving points simultaneously, which corresponds to about the same computing time as for onereceiving point using the reciprocity principle method. Because the interface was representedusing many small segments rather than using square or rectangular cell boundaries, the traveltimescalculated are more accurate than those using the reciprocity principle method.536.4 Effects on the ComputationThe effects of cell shape, cell size, and the magnitude of the velocity field were tested byusing a few simple models. The effects of the grid size were tested using model 2, a velocityfield with a constant vertical gradient. The traveltimes computed by the Vidale method andthe NLTI method were compared with the true times. As expected, the errors decreased forboth methods when a denser grid was used. A constant velocity model was also used to testthe effects of using a rectangular grid. Square cells were substituted by rectangular cells byincreasing the grid sampling rate in the vertical direction. For some of these rectangular cells,the traveltimes of two grid points at a shorter cell boundary were used to calculate the timesof the other two grid points. Since the wave fronts for shorter cell boundaries were closer toplane waves, traveltime accuracies were improved in these regions. Tests on the effects of themagnitude of the velocity field showed that the computation errors of the NLTI method weresmaller in a high velocity field as compared to a low velocity model.6.5 Application in Traveltime TomographyThe FD methods for first arrival traveltime computation which exist in the literature arethe Vidale method and the LTI method. No interfaces were modelled in these methods. Thealgorithm developed in this thesis extended the first arrival traveltime computation to a layeredmodel with an arbitrary velocity distribution in each layer. This algorithm can also be usedto compute reflection traveltimes. It is considered to be faster and more accurate than thereciprocity principle method. Therefore, it provides an efficient tool in solving the forwardproblems in transmission and reflection traveltime tomography.54REFERENCESAldridge, D. F., and Oldenburg, D. W., 1992, Two dimensional tomographic inversion with finitedifference traveltimes, Submitted to J. of Seis.Aldridge, D. F., 1992, Comment on "Conversion points and traveltimes of converted waves inparallel dipping layers" by Gisa Tessmer and Alfred Behle, Geophysical Prospecting, 40,101-103.Asakawa, E., and Kawanaka, T., 1991, Seismic ray tracing using linear traveltime interpolation,unpublished.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, H. W., 1985, Tomographic determination of velocity anddepth in laterally varying media, Geophysics, 50, 903-923.Bording, R. P., Gersztenkorn, A., Lines, L. R., Scales, J. A. and Treitel, S., 1987, Applicationof seismic travel-time tomography, Geophys. J. R. astr. Soc., 90, 285-303.Dix, C. H., 1981, Seismic Prospecting for Oil, 2nd edition, Intern. human resour. devel. corp.,p125-126. pterven, V., Molotkov, L A. and Pgen6ik, I., 1977, Ray method in seismology, UniverzitaKarlova, Praha.Itoh, K., Saito, T. and Yamada, H., 1983, A method for approximate computation of first arrivaltraveltimes in seismic refraction modeling, Proc. 1983 Spring Meeting, SEG Japan, 11-12,(in Japanese).Mandal, B., 1992, Ray Tracing and Travel Time Using Triangular Grid, 54th Meeting andTechnical Exhibition, European Association of Exploration Geophysicists, 498-499.Matsuoka, T., Asakawa, E., and Kawanaka, T., 1990, Forward modeling for ray tomograghy,Geotomography, vl, The First SEGJ International Symposium on Geotomography, 185-195.Matsuoka, T., and Ezaka, T., 1990, Ray tracing on reciprocity, SEG 1990 Meeting, 1028-1031.55Matsuoka, T., and Ezaka, T., 1992, Ray tracing using reciprocity, Geophysics, 57, 326-333.Moser, T. J., 1989, Efficient seismic ray tracing using graph theory, SEG 1990 Meeting, Expandedabstracts, 1106-1108.Moser, T. J., 1991, Shortest path calculation of seismic rays, Geophysics, 56, 59-67.Podvin P., and Lecomte, I., 1991, Finite difference computation of traveltimes in very contrastedvelocity models: a massively parallel approach and its associated tools, Geophys. J. Int.,105, 271-284.Qin, F., Olsen, K. B., Luo, Y., and Schuster, T., 1990, Solution of the eikonal equation by afinite-difference method, SEG 1990 Meeting, 1004-1007.Qin, F., Luo, Y., Olsen, K. B., Cal, W., and Schuster, T., 1992, Finite difference solution of theeikonal equation along the expanding wavefronts, Geophysics, 57, 478-487.Saito, H. and Ohtomo, H., 1986, An attempt to reduce the occurrence of the artifacts in seismicray tomograph, 56th annual international SEG Meeting, Houston.Saito, H., 1989, Traveltimes and raypaths of first arrival seismic waves: computation methodbased on Huygens' principle, SEG 1989 Meeting.Saito, H., 1990, 3—D ray-tracing method based on Huygens' principle, SEG 1990 Meeting.Schneider, W. A., Ranzinger, K. A., Balch, A. H., and Kruse C., 1992, A dynamic programmingapproach to first arrival traveltime computation in media with arbitrarily distributedvelocities, Geophysics, 57, 39-50.Telford, W. M., Geldart, L. P., Sheriff, R. E. and Keys, D. A., 1988, Cambridge UniversityPress, 271-276.Thurber, C. H. and Ellsworth, W. L., 1980, Rapid solution of ray tracing problems inheterogeneous media, Bull. Seis. Soc. Am., 70, 1137-1148.Trier J. V., and Symes W. W., 1991, Upwind finite-difference calculation of traveltimes,Geophysics, 56, 812-821.Vidale, J., 1988, Finite-difference calculation of traveltimes, Bull. Seis. Soc. Am., 78, 2062-2076.56White, D. J., 1989, Two-dimensional seismic refraction tomography, Geophysical Journal, 97,223-245.57Appendix APossible Interface-cell CombinationsAn interface segment can cross a cell in twenty possible ways which are illustrated inthe diagram below. During the discretization of the model the intercept point coordinatesare computed. The irregular cell types, the number inside the parentheses, and the irregularcell parameters, e, el and e2 in the diagram, are memorized by the program and are used todecompose the irregular cell into smaller rectangular or triangular cells.(1) (2) (3) (4)(5) (6) e 2 (8)e 2(9) (10) (11) (12)e 2 ei e2e 2(13)e 2(14) (15) (16)e 2/e(17) (18) (19)^(20)58Appendix BTraveltime Computation for a Point in a Vertically Increasing Velocity ModelSuppose that, for a vertical gradient velocity medium, v(x, z) = vo kz, where vo is thevelocity at zero depth and k is the velocity gradient, s is the shot point and also the origin ofthe coordinates, and point A, with the coordinate (x, z), is the point at which the traveltime isto be computed. Then, the following can be shown to hold (Telford et al., 1988, 271-276):1. the ray paths are the arcs with the centers of the circles located at line z = —vo/k;2. a wave front is an arc with center of the circle at the vertical axis and radius r 2 = CS. CS';3. the maximum depth a ray can reach is hinaz = V k(COSh kt — 1), the offset of themaximum depth point r max = vo/ksinh kt and the time reach the maximum point tmax1iksinh -1 (krmax/vo).Using the above, the traveltime for any point A(x, z) can be computed. Suppose the wavefront circle passing point A is centered at point C, and BC is a horizontal line segment withits length equal to the wave front radius r. There must exist one ray circle with B as the59maximum depth point, i.e., the wave front circle and ray circle encounter at point B in rightangle — one in vertical and one in horizontal. We can relate r to the coordinates of point Aby r2 = (z — h)2 + x2. Also, r2 h(h 2v0/k). Eliminating r, we obtain the relationshipbetween h and (x, z) as h kl2(x2 z2)I(vo kz). Since h = volk(coshkt — 1), we obtainthe traveltime for point A(x,z) ast (x, z) —1 sinh- 1 (21)
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Geometrical formulation and traveltime computation...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Geometrical formulation and traveltime computation in heterogeneous layered media Li, Xingong 1992
pdf
Page Metadata
Item Metadata
Title | Geometrical formulation and traveltime computation in heterogeneous layered media |
Creator |
Li, Xingong |
Date Issued | 1992 |
Description | This thesis deals with the subject of seismic traveltime computation. The velocity field is represented by an arbitrary distribution of rectangular or square cells. The cells with an interface crossed are represented by two non-rectangular cells. Each cell is assigned a constant velocity. A geometrical approach is presented to construct the Vidale formula (Vidale, 1988) and the linear traveltime interpolation, LTI, formulae (Matsuoka et al., 1990; Asakawa and Kawanaka, 1991; Podvin and Lecomte, 1991; Schneider et al., 1992). The geometrical construction allows an easy visualization of the principles which underlie both the Vidale and the LTI methods. From this geometrical point of view, subdividing the cell boundary into smaller segments (the subdivided LTI, or the SLTI method) improves the approximation of the plane wave assumption and hence, improves the resulting accuracy. A new approach is developed to calculate the traveltimes for a heterogeneous model with interfaces. The algorithm involves three steps: 1. The discretization of the velocity field using square or rectangular cells. For cells on an interface, each cell is divided into two secondary non-rectangular cells. The two non-rectangular cells are assigned different velocities, so that even an interface with a large contrast in velocities can be modelled. 2. The calculation of traveltimes at the interface locations. 3. The computation of the reflection traveltimes and the transmission traveltimes by expanding the wave fronts from the interface points upwards and downwards from the interface. The principles used in this approach are Fermat's principle and Huygens' principle. The known grid points or interface points are treated as secondary sources and the traveltimes at the points within a cell are calculated and revised by finding the shortest traveltimes from the known points. Because of the special consideration given to the interface cells, the new approach is found to be more accurate than the Vidale and the LTI approaches. With respect to reflection traveltime computation, the new algorithm calculates the reflection traveltimes for all the receiving points at the same time and is, consequently, faster than the reciprocity principle method (Matsuoka and Ezaka, 1991; Podvin and Lecomte 1991). Various examples using the algorithm are presented. Model experiments show that the computation errors are smaller in a high velocity field than in a low velocity field. The rectangular grid is tested by increasing the vertical sampling rate of the square grid by a factor of two. The results show that computational errors decrease in the horizontal direction and increase in the vertical direction. |
Extent | 2535741 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-07-29 |
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.0052959 |
URI | http://hdl.handle.net/2429/1186 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1993-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1993_spring_li_xingong.pdf [ 2.42MB ]
- Metadata
- JSON: 831-1.0052959.json
- JSON-LD: 831-1.0052959-ld.json
- RDF/XML (Pretty): 831-1.0052959-rdf.xml
- RDF/JSON: 831-1.0052959-rdf.json
- Turtle: 831-1.0052959-turtle.txt
- N-Triples: 831-1.0052959-rdf-ntriples.txt
- Original Record: 831-1.0052959-source.json
- Full Text
- 831-1.0052959-fulltext.txt
- Citation
- 831-1.0052959.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-0052959/manifest