Published on: **Mar 3, 2016**

- 1. Mathematical representation models and applications on seismic tomography M. Arvanitis and B. D. Al-Anazi REVIEW New mathematical techniques have contributed substantially to the improvement of the geophysical prospecting methods, like traveltime seismic tomography. Thanks to these new techniques, the time to solve an inverse problem has been reduced dramatically making seismic tomography applicable to a great number of problems even in three dimensions. New raytracing and wavefront techniques provide a more flexible parameterization. Advancement from the least squares technique to today’s back-projection method’s, for example, has improved tomographic methods. Key words: tomography, raytracing, wave front, grid points, travel time INTRODUCTION What does tomography mean as a word? Its origin is found in the Greek word “tomo”, which means slice. Therein lies the basic idea: if we take many 2-D slices, ac- cording to the central slice theorem28, we can reconstruct the whole 3-D image of an object. Thanks to the same theorem, we can easily construct 2-D sections from 1-D lines, which can be measured in experiments. Seismic travel time tomography can be defined as the reconstruction of the Earth’s velocity model, using the seismic waves travel time deviations from a reference ve- locity model, better known as starting or background model. A starting model is an initial guess, an estimate of a velocity model. Seismic tomography as we know it to- day originated in 1974 as “3-D inversion”. At first, seis- mologists were very skeptical about the new method and its results. The whole attitude changed dramatically in the mid-1980’s, when iterative methods were intro- duced, to facilitate the calculation of large and sparse matrixes that occurred from the seismological data.6,23 Believability in the method was linked to the first global tomographic results9,37,11,which correlated satisfactory with the geoid. As credibility of the method grew seismol- ogists renamed “3-D inversion” to what is known today as seismic tomography. We can classify seismic tomography in two main categories; travel time and amplitude tomography. In this paper we will focus only on travel time tomography. Regarding the nature of the seismic waves, travel time tomography can be divided into refraction, reflection and diffraction tomography. Referring to the source, whether it is a natural earthquake or a shot, we carve up tomogra- phy into passive and active tomography, respectively. Seismic tomography is an imaging technique that uses seismic waves generated by earthquakes and explosions to create computer-generated, three-dimensional images of Earth's interior. This is how seismologists infer the different layers in the Earth. How is this done? The time it takes for a seismic wave to arrive at a seismic station from an earthquake can be used to calculate the speed along the wave's ray path. By using first arrival times of P waves recorded by seismic stations all over the world, scientists are able to define slower or faster regions deep in the Earth The simplest case of seismic tomography is to estimate P-wave velocity. Several methods have been developed for this purpose, e.g., refraction traveltime tomography, fi- nite-frequency traveltime tomography, reflection traveltime tomography, waveform tomography.36 To obtain a higher-resolution image one has to abandon the infinite-frequency approximations of ray theory that are applicable to the time of the wave 'onset' and instead measure travel times (or amplitudes) over a time window of some length using cross-correlation. Finite-frequency tomography takes the effects of wave diffraction into ac- count, which makes the imaging of smaller objects or anomalies possible.24 The raypaths are replaced by volumetric sensitivity kernels, often named 'banana-doughnut' kernels in global tomography, because their shape may resemble a banana, whereas their cross-section looks like a doughnut, with, at least for direct P and S waves, zero sensitivity of the travel time on the geometrical ray path. In finite-frequency tomography, travel time and amplitude anomalies are frequency-dependent, which leads to an increase in resolution. To exploit the information in a seismogram to the full- est, one uses waveform tomography. In this case, the seismograms are the observed data. In seismic explora- tion, the forward model is usually governed by the acous- tic wave equation. This is an approximation to the elastic wave propagation.36 Elastic waveform tomography is much more difficult than acoustic waveform tomogra- phy. The acoustic wave equation is numerically solved by some numerical schemes such as finite-difference and fi- nite-element methods. In global tomography the inverse problem for elastic waves can be handled by adjoint methods. PARAMETERIZATION Lets consider two closely spaced points in the medium; the inverse of the local wavespeed associated with these NAFTA 60 (9) 495-498 (2009) 495
- 2. points is the slowness. There are three kinds of slowness models: homogeneous and heterogeneous cells (2-D) or blocks (3-D) of constant slowness values and rectangular grids with slowness values assigned to the grid points with different interpolation schemas to specify the values between the grid points. In general, the use of cells is the most common parameterization but it is facing difficul- ties, since the sharp boundaries between the cells cannot be resolved. As for the grid parameterization, a fine regu- lar and irregular grid parameterization exist. The former parameterization is the purely tomographic approach while the latter one is closer to forward modelling. Regu- lar grid has the advantage of simplicity but it can cause over-parameterization when high resolution is required. Recent studies are focused on irregular grid using Delaunay triangles or Veronoi polygons to avoid such problems.4,41 The travel time for a ray is: ( )T s V r s= ò 1/ d (1) where V(r) is the unknown velocity for the ray-path S. We want to determine V(r) from N travel time measure- ments. Let To be the travel time for the starting model: ( )T s V r s0 0 01= ò / d (2) Whether we are not sure about the estimate of the start- ing model we can use a tau-p method from picked arrival times to form a reliable starting velocity model.2 Using Fermat.s principle, we can ignore the true ray-path and use the ray-path of the starting model instead of it. The delay time is: ( ) ( ) ( )dT T T s V r s s V r s s V V s= - = - » - »ò òò0 0 0 0 01 1 1 1/ / / /d d d ( ) ( )( )» - -ò s V r V r s0 0 2 d d (3) where: ( ) ( ) ( )dV r V r V r= - 0 (3a) Equation (3) comprises a linear system of equations, which can be changed in such a way to become more suit- able for computer processing. We parameterize the me- dium with I interpolation functions hi, which is the basis of the subspace of the Hilbert space of all possible mod- els V(r): ( ) ( )d gV r h rk k k= å (4) where k spans the integers from 1 to I, and function gk is the weight of the function hk. Considering a cell- parameterization we have: hi (r)=1, if r in cell I and hi (r)=0, anywhere else (5) Equation (3) can now be defined as: ( ) ( )[ ]d g gT s V r h r s Ak k k ik k k= - =òå å0 0 2 1/ d (6) where: ( ) ( )[ ]A s V r h r sik k= -ò 0 0 21/ d (6a) Equation (6) can be formulated for each shot in a ma- trix form and in terms of slowness as: Ms t= (7) Where s is the slowness vector, t the time vector and M is the matrix of lij, where lij is the length of the i-th ray-path through j-th cell. SOLUTION Foremost, we have to calculate the matrix elements Aik, which implies the finding of the ray-path. Two methods are commonly used in seismic tomography to find the ray path: ray-tracing and wavefront methods. The two-point ray-tracing finds ray-paths along which seismic energy is propagating and calculates the travel time. For a layered media, rays are traced by solving the differential equa- tions with a Runge-Kutta predictor-corrector scheme. To define in a better way the ray geometry and the slowness, we mainly have two methods: shooting and bending, in- spired both from ray-tracing. The former is based on continuous iterations until the end of a ray to meet a limit condition or by interpolating between close rays using hermite cubic interpolation.7 The latter uses a parameterization for a ray-path by the support points Vi of a third order B-spline. The location of the ray is a func- tion of the four nearest points: ( )Q u b V b V b V b Vi i i i= + + +- - - - +2 2 1 1 0 1 1 (8) Where b1 depends on u, 0 £ u £ 1, and is known for the different values of u. A conjugate gradient method33 is needed to find which of the support points V1 minimizes the time given by (1). Èervený5 uses the quadratic slow- ness, 1/V2, instead of slowness since it offers the simplest analytical solution in inhomogeneous medium. Latter techniques pored over the drawbacks of ray-tracing in- cluding the works of Zelt and Ellis49 who invented a ray-tracing technique with a trapezoidal parameteri- zation providing a rapid traveltime calculation and of Sethian and Popovici34 who presented the fast marching technique that can model turning rays, but with unsatis- factory accuracy. The wavefront (surfaces of equal travel time) methods are alternative techniques to ray-tracing. This method determines minimum ray-paths and traveltimes by ex- panding a wavefront in the whole model. Recent ad- vances in the ray-tracing methods have mainly been focused on wavefront methods rather than ray-tracing, for two basic reasons: ray-tracing is valid only for smooth velocity structures and it is significantly slower than any wavefront method. Vidale43,42 modified the wavefront method and the eikonal equation solver introducing a fi- nite difference procedure, to propagate traveltimes through a uniformly sampled grid. The eikonal solver finds the wavefront that forms concentric shells about the source and conducts the ray-paths from their shape. A defect of Vidalia’s method27 is that it fails when velocity contrasts are of the order of u u2 1 2/ > (9) (Hole et al.14 modified Vidalia’s algorithm to a more rapid algorithm using variable grid spacing. SPR, for 496 NAFTA 60 (9) 495-498 (2009) M. ARVANITIS AND B. D. AL-ANAZI MATHEMATICAL REPRESENTATION MODELS AND APPLICATIONS...
- 3. Shortest Path Ray-tracing, developed primarily by Saito32,31 and Moser20, expands the wavefront in the en- tire velocity network being in that way more stable to ve- locity contrasts. Zhang, et al.53,54 developed a SPR method with a graphic template that allowed only straight rays within a constant cell, calculating rapidly and with good accuracy ray-paths and traveltimes on a large number of grids. In the wavefront construction44,45 the entire wavefront is represented by a triangular mesh providing good accuracy but implying time-consuming calculations. GRT, for grid ray-tracing55, combines the advantages of both wavefront construction and fast marching, by tracing rays within a local grid. GRT is about 8 times faster than any wavefront method. The main shortcoming of the wavefront methods is that they used to calculate only first arrivals. Moser20, Hole and Zelt13 and Zhang and Toksöz54 modified the method to calculate also latter arrivals. The solution of the equation (7) is challenging. With real data there always exists possibility that no ray crosses a cell (an ill-posed system). The other problem is the amount of data that we use in iterative solution tech- niques or inversion. Early tomographic methods used exact least squares techniques to solve the resulting sys- tem of equations, but model limitations (due to the re- strained computer capabilities) on the order of 1 000 unknowns were coming up. This number is not that large or as much of a concern today, e.g. the tomographic anal- ysis of Zelt and Barton47 involves more than 50 000 travel times. To surpass these obstacles, iterative matrix solvers were introduced. McMehan18 and Neumann- Denzau and Behrens22, improved an older iterative method of Kaczmarc and renamed it to ART, for Alge- braic Reconstruction Technique. This method was badly conditioned and extremely slow for seismic tomography. Gilbert10 invented a more efficient method, named SIRT, for Simultaneous Iterative Reconstruction Technique. Both ART and SIRT are mainly applicable when pixels or voxels (the three dimensional analog of pixels) are used as the basis function (5). Both of these techniques use backprojection in an iterative manner to solve the system of equations. Backprojection is an iterative process to es- timate the average slowness. Instead of backprojecting travel time residuals along ray-paths, other formulations backproject phase residuals along wavepaths35 which take into account the finite frequency effects in travel time data. In order to reconstruct velocities and interfaces we solve the regularized inverse problem. The trial-and-error forward modelling is a time-consuming and laborious process which fails to provide an estimate of parameter uncertainty and resolution, as inversion does. In general, the non-linear system of equations (7) is being linearized, considering that the velocity structure is divided into a reference or starting model, that is as- sumed to be known, and an unknown perturbation, which is considered as very small. When the model is non-linearized, the solution will be independent of the model parameterization, so we apply the Tikhonov method40 in order to reconstruct the model with a Laplacian operator. Broadly speaking, for a regularised problem we can invert various types of data, e.g. reflec- tion, refraction, crosswell, both refraction and reflection data, known as joint inversion, for better performance and various types of model parameters (e.g. slowness, re- flector geometries). The aim of the regularized inversion is to minimize a tradeoff parameter concerning the data misfit. In early tomographic problems fitting data was the major care, but recent studies prove that we can fit data to any small misfit magnitude according to the con- straints of the model parameters, although the solution may not be physically consistent.38 Given that a crucial demerit of inversion is its nonuniqueness30, which leads to multiple solutions of the problem, the main concern is how to obtain a stable and unique solution that doesn't provide unneeded structures. Fitting travel times with a least squares criterion can't always provide the best solu- tion, although the vast majority of the tomographic meth- ods implement some variant of the principal least-square method, by selecting a model that minimizes a certain measure of travel time error, e.g. the damped least-squares technique17, the conjugate gradient meth- ods, the biconjugate gradient algorithm25,26 etc. Another approach54 is to regularize average slowness (traveltime divided by ray-length) and apparent slowness (traveltime derivative with respect to surface distance) than traveltimes or to apply smoothness constraints or deriv- ative operators to find the simplest structure that fits the data under a given tolerance. A generally accepted opinion is that 2-D seismic inver- sion can give an incorrect picture of the subsurface21 so 3-D inversions are needed in order to put up a better im- age. In the last two decades, three dimensional seismic methods have been developed but early studies, (Thurber39, Kanasewich and Chiu15,) faced many prob- lems due to limitations in computer resources and data coverage. Zelt47 (1994) inverted simultaneously refrac- tion and reflection travel time data in order to provide a three dimensional starting model. The recent develop- ments in 3D refraction methods48,3,8,12 have partially solved many problems of the past, but there is still a long way to go, e.g. difficulties still occur when large velocities contrasts exist, wavespeeds cannot be resolved reliably in the main refractor.16 CONCLUSIONS With the help of new mathematical techniques, seismic tomography has been evolved to a widely used technique covering a broad range of applications from global to- mography to near surface geophysics. The basis of a tomographic problem is the inversion of a matrix. The re- quested precision as well as today’s technical needs de- mand us to solve systems of the order of more than, e.g. 10 000 unknowns. The inversion of such matrices can- not be reliably handled with conventional techniques, so new inversion techniques like, e.g. SIRT, have been intro- duced. Today, to deal with similar problems, modern backprojection techniques are used. Lately, many efforts have been focused on more advanced methods, like the inversion using genetic algorithms. MATHEMATICAL REPRESENTATION MODELS AND APPLICATIONS... M. ARVANITIS AND B. D. AL-ANAZI NAFTA 60 (9) 495-498 (2009) 497
- 4. REFERENCES 1. Aki K., Christofferson, A., Husebye, E.S., et al., Three-dimensional seismic velocity anomalies in the crust and upper mantle under the U.S.G.S California seismic array, EOS, Trans. Am. Geophys. Union, 56, 1145, 1974. 2. Barker, N., Barton, P. and Nicola-Carena, E., Rapid automatic velocity images in depth and time from long-offset two-ship datasets, LITHOS science report 2, 35- 42, 2000. 3. Bennet, G., 3D seismic refraction for deep exploration targets, The Leading Edge, 18, 186-191,1999. 4. Bohm, G., Galuppo, P., and Vesnaver, A.A., 3D adaptive tomography using Delaunay triangles and Voronoi polygons, Geophysical Prospecting, 48, 723-744, 2000. 5. Èervený, V., Ray-tracing algorithms in three-dimensional laterally varying lay- ered structures, ed. D.Reidel, Norwell, Mass, 1987. 6. Clayton, R.W. and Comer, R.P., A tomographic analysis of mantle heterogeneities from body wave travel time, EOS, Trans. Am. Geophys. Un- ion, 64, 776, 1983. 7. Comer, R.P., Rapid seismic ray-tracing in a spherically symmetric Earth via in- terpolation of rays, Bull. Seism. Soc. Am., 74, 479-92, 1984. 8. Deen, T.J., Gohl, K., Leslie, C., Papp, E., Wake-Dyster, K., Seismic refraction inversion of a paleochannel system in the Lachlan Fold Belt, Central New South Wales, Exploration Geophysics, 31, 389-393, 2000. 9. Dziewonski, A.M. and Anderson, D.L., Seismic tomography of the Earth.s in- terior, Am. Sci., 721, 483-494, 1984. 10. Gilbert, P.F.C., Iterative methods for three dimensional reconstruction from projections, J. Theor., Biol., 36, 105-17, 1972. 11. Hager, B.H., Clayton, R.W., Richards, M.A., et al., Lower mantle heterogene- ity, dynamic topography and the geoid, Nature, 313, 541-545, 1985. 12. Hobro, J., Seismic traveltime tomography in three dimensions: synthetic ex- periments and results obtained using Jive3D, LITHOS science report, 2, 11-20, 2000. 13. Hole, J.A., and Zelt, B.C., 3D finite-difference reflection travel times, Geo- physical Journal International, 121, 427-434, 1995. 14. Hole, J.A., Clowes, R.M., Ellis, R.M., Interface Inversion using broadside seismic refraction data and three dimensional travel time calculations, Jour- nal of Geophysical Research, 97, 3417-3429, 1992. 15. Kanasewich, E.R., Chiu, K.L., Least squares inversion of spatial seismic re- fraction data, Bull. Seismol. Soc. Am., 75, 865-880, 1985. 16. Lanz, E., Maurer, H., and Green, A.G., Refraction tomography over a buried waste disposal site, Geophysics, 63, 1414-1433, 1998. 17. Lutter, W.J., Nowack, R.L., Inversion for crustal structure using reflections from the PASSCAL Ouachita experiment, J.Geophys.Res., 95, 4633-4646, 1990. 18. McMehan, G.A., Seismic tomography in boreholes, Geophys.J. R. Astron. Soc., 74, 601-612, 1983. 19. Menke, W., Geophysical data analysis: discrete inverse theory, Academic Press, 1984. Moser, T.J., Efficient seismic ray-tracing using graph theory, 59th Annual International Meeting, SEG, 1989. 20. Moser, T.J., Shortest path calculation of seismic rays, Geophysics, 56, 59-67, 1991. 21. Nestvold, E.O., 3-D seismic: is the promise fulfilled?, The Leading Edge, 11, 12- 19, 1992. 22. Neumann-Denzau, G. and Behrens, J., Inversion of seismic data using tomographic reconstruction techniques for investigation of laterally inhomogeneous media, Geophys. J. R. Astron., Soc., 79, 305-316, 1984. 23. Nolet, G., Solving or resolving inadequate and noisy tomographic systems, J. Comput. Phys., 61, 463-482, 1985. 24. Nolet, G., "A Breviary of Seismic Tomography", Cambridge University Press, 2008 25. Paige, C.C., and Saunders, M.A., LSQR: An algorithm for sparse linear equa- tions and sparse least squares, ACM Trans. Math. Software, 8, 43-71, 1982. 26. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, Numerical recipes in Fortran: the art of scientific Computing, Cambridge University Press, 1992. 27. Qin, F., Luo, Y., Olsen, K.B., Cal, W., Schuster, G.T., Finite-difference solu- tion of the eikonal equation along expanding wavefronts, Geophysics, 57, 478-487, 1992. 28. Radon, J., Über die Bestimmung von Fonktionnen durch ihre Integralwerte langs gewisser Mannigfaltigkeiten, Ber.Verh. Saechs. Akad. Wiss, Liepzig. Math. Phys. Kl., 69, 262-277, 1917. 29. Ramillien, G., Genetic algorithms for geophysical parameter inversion from altimeter data, Geophysics Journal International, 147, 393-402, 2001. 30. Ramm, A.G., Examples of non-uniqueness for an inverse problem of Geo- physics, Appl. Math. Lett, 8, 87-90, 1995. 31. Saito, H., 3-D ray-tracing method based on Huygens. principle, 60th Annual International Meeting, SEG, 1990. 32. Saito, H., Traveltimes and ray-paths of first arrival seismic waves: computation method based on Huygens. principle, 59th Annual International Meeting, SEG, 1989. 33. Scales, J., Tomographic inversion via the conjugate gradient method, Geo- physics, 52, 179-185, 1987. 34. Sethian, J.A., and Popovici, A.M., 3-D traveltime computation using the fast marching method, Geophysics, 64, 516-523, 1999. 35. Sheng, J. and Schuster, G., Finite-Frequency resolution limits of wavepath traveltime tomography for smoothly varying velocity models, 2002, SEG An- nual Meeting 2002. 36. Stewart, R. R., Exploration Seismic Tomography: Fundamentals, Society of Exploration Geophysicists, 1991 37. Tanimoto, T. and Anderson, D.L., Mapping convection in the mantle, Geophys. Res. Lett., 11, 287-290, 1984. 38. Thompson, D., Nonlinear waveform tomography, Ph.D Thesis, MIT, Cam- bridge, 1993. 39. Thurber, C.H., Earthquake location and three-dimensional crustal structure in the Coyote Lake Area, central California, J. Geophys. Res, 88, 8226-8236, 1983. 40. Tikhonov, A.N. and Arsenin, V.Y., Solutions of ill-posed problems, Winston, 1977. 41. Trinks, I., Singh, S., Chapman, C., Barton, P., Bosch, M., Traveltime tomogra- phy using irregular parameterised grids, Journal of Sub-basalt conference ab- stracts, 7, 193, 2002. 42. Vidale, J.E., Finite-difference calculation of travel times in three dimensions, Geophysics, 55, 521-526, 1990. 43. Vidale, J.E., Finite-difference calculation of travel times, Bulletin of the Seis- mological Society of America, 78, 2062-2076, 1988. 44. Vinje, V., Iversen, E., Gjoystdal, H., Traveltime and amplitude estimation us- ing wavefront construction, Geophysics, 58, 1157-1166, 1993. 45. Vinje, V., Iversen, E., Gjoystdal, H., and Astebrl, K., Estimation of multivalued arrivals in 3D models using wavefront construction, Geophys. Prospecting,, submitted 1996, in press. 46. Zelt, C.A and Zelt, B.C., Study of out-of-plane effects in the inversion of re- fraction/wide angle reflection traveltimes, Tectonophysics, 286, 209-221, 1998. 47. Zelt, C.A. and Barton, P.J, 3D seismic refraction tomography: a comparison of two methods applied to data from the Faroe Basin, Journal of Geophysical re- search, 103, 7187-7210, 1998. 48. Zelt, C.A., 3-D velocity structure from simultaneous traveltime inversion of in-line seismic data along intersecting profiles, Geophysical Journal Interna- tional, 118, 795-801, 1994. 49. Zelt, C.A., and Ellis, R.M., Practical and efficient ray-tracing in two dimen- sional media for rapid traveltime and amplitude forward modelling, Canadian Journal of Exploration Geophysics, 24, 16-31, 1988. 50. Zelt, C.A., Lateral velocity resolution from 3D seismic refraction data, Geo- physical Journal International, 135, 1101-1112, 1998. 51. Zelt, C.A., Modelling strategies and model assessment for wide-angle seismic travel time data, Geophysical Journal International, 139, 183-204, 1999. 52. Zelt, C.A., Smith, R.B., Seismic travel time inversion for 2D crustal velocity structure, Geophysics Journal International, 108, 16-34, 1992. 53. Zhang, J., Rodi, W., Mackie, R.L., Shi, W., Regularization in 3-D DC resistiv- ity tomography, presented at SAGEEP, 1996. 54. Zhang., J., ten Brink, U.S., Toksöz, M.N., Nonlinear refraction and reflection travel time tomography, Journal of Geophysical Research, 103, 29743-29757, 1998. 55. Zhu, T., and Cheadle, S., A grid raytracing method for near-surface traveltime, 69th Annual Meeting SEG, Expanded Abstracts, 1759-1763, 1999. v Authors: Michael Arvanitis, Geomorph Instruments, Greece, mike@geomorph.gr Bandar Duraya Al-Anazi, King Abdulaziz City for Science & Technol- ogy,Saudi Arabia, bandar.alanazi@gmail.com 498 NAFTA 60 (9) 495-498 (2009) M. ARVANITIS AND B. D. AL-ANAZI MATHEMATICAL REPRESENTATION MODELS AND APPLICATIONS...