Published on: Mar 3, 2016
Transcripts - nafta mathemat
Mathematical representation models and
applications on seismic tomography
M. Arvanitis and B. D. Al-Anazi
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
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
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
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
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
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
d d (3)
( ) ( ) ( )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-
( ) ( )d gV r h rk
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= - =òå å0 0
1/ d (6)
( ) ( )[ ]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.
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-
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...
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
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
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
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,
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,
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,
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,
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,
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,
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,
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-
39. Thurber, C.H., Earthquake location and three-dimensional crustal structure
in the Coyote Lake Area, central California, J. Geophys. Res, 88, 8226-8236,
40. Tikhonov, A.N. and Arsenin, V.Y., Solutions of ill-posed problems, Winston,
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,
55. Zhu, T., and Cheadle, S., A grid raytracing method for near-surface
traveltime, 69th Annual Meeting SEG, Expanded Abstracts, 1759-1763, 1999.
Michael Arvanitis, Geomorph Instruments, Greece, email@example.com
Bandar Duraya Al-Anazi, King Abdulaziz City for Science & Technol-
ogy,Saudi Arabia, firstname.lastname@example.org
498 NAFTA 60 (9) 495-498 (2009)
M. ARVANITIS AND B. D. AL-ANAZI MATHEMATICAL REPRESENTATION MODELS AND APPLICATIONS...