Determination and sensitivity analysis of the seismic velocity of a shallow layer from refraction traveltimes measures

In this paper, we are interested in determining the seismic velocity of a shallow under-ground layer from refraction traveltimes measures. We present a study case taken from an experimental seismic survey. The study case is a wide-angle seismic inversion using experimental traveltimes measures and based on ray tracing technique and genetic algorithms. The hypothesis on the velocity distribution, coming from the seismic experiment, makes the computation of some seismic rays expensive in time. We propose to reduce the computations time by introducing a formulation of the inverse problem that avoids such costly rays, hence the inversion becomes feasible. Also we present a sensitivity analysis based on a singular value decomposition of the jacobian of the traveltimes with respect to velocity. We give the relationship between the traveltimes measure errors and the velocity estimation error. We discuss the advantages of this method over the classical one based on the resolution matrix.


INTRODUCTION
Geological investigations are usually based on seismic surveys that allow obtaining images of under-surface sublayers and interfaces between these layers.One of the aims of these surveys is to carry out an estimation of the depth of the bedrock and the mechanical properties of the rocks.This domain of geophysics is referred to "seismic inversion" as it deals with seismic traveltimes measures to recover the seismic velocity in the underground.Knowing its industrial applications, this tool has become in recent years very powerful to reach high resolution images for complex geological structure.
In this paper, we are interested in carring out a seismic refraction traveltimes inversion to determine the seismic wave's velocity in a shallow underground sublayer.The inverse problem is about seismic refraction in a two dimensional domain composed of two layers (Ω 1 , Ω 2 ) and separated by an interface (Γ).Let (v 1 , v 2 , γ) be the velocity in (Ω 1 ), the velocity in (Ω 2 ) and the interface between them respectively.We call the vector of these three functions the continuous velocity distribution.Seismic traveltimes T are related to the wave's velocity through the Eikonal equation: (1) where υ(x) is the wave's velocity at point x of the domain.Highlights of this model are presented in section 2.
A seismic traveltime inversion seeks to find the distribution of the wave's velocity (υ 1 , υ 2 , γ) in the underground from an array d of N T experimental traveltime measures corresponding to some positions of the seismic source s i and the receiver r i , i = 1 . . .N T .The choice of the model space and the data space, to which belong the velocity distribution and the traveltimes measures, will be presented in section (3).The inverse problem is usually formulated as a nonlinear least squares optimization problem with a cost function (d -BT (υ)) t C d -1 (d -BT (υ)) being the misfit function between the vector of experimental data d and simulated traveltimes vector BT (υ) for a given seismic velocity υ, B is the operator that gives an array of values of T (υ) for the pairs (s i , r i ) and C d is the data covariance matrix (see [1], [2] and [3]).The nonlinear least squares problem is solved using stochastic ("genetic algorithms") and/or deterministic methods.The literature associated with these two classes of methods is huge and no attempt to review it shall be made in this paper.
In both cases it is assumed that, for any velocity distribution (υ 1 , υ 2 , γ), the Eikonal solver must provide a vector BT (υ) of N T elements to be able to calculate the misfit function.In many applications this condition is trivial.However, when the depth of the interface γ is very smaller then the distance between the source and the receiver (shallow layer), the computation of all the components of the vector BT (υ) may be computationally not feasible.The reasons of this are elucidated in paragraph (2.2).To overcome this difficulty we propose in section 5 a formulation of the inverse problem that avoids computing expensive traveltimes, thus allowing the inversion become feasible.
The sensitivity analysis of a mathematical model, describing a physical system, is the study of how errors in the output data (traveltimes measures) are affected by the perturbation of input parameters (seismic velocity).The process of producing a solution to an inverse problem is not complete without a sensitivity analysis of the solution with respect to the measure errors.Thus in addition, one must study the robustness or the quality of the solution, and not simply producing a single solution that minimizes an objective function.A common quantitative approach to study the sensitivity of the solution is the computation of the diagonal elements of the resolution matrix.As explained in [4], [5] and [6], the expression of the resolution matrix is derived by assuming that the traveltime measures are noise free.This is why this commonly used method does not provide a direct relation between the estimation error of the velocity distribution and the measure errors.
In this paper, we consider a class of sensitivity analysis methods called optimal vectors method.It is based on computing the singular value decomposition (SVD) of the jacobian of the output data with respect to the input parameters (see [7] and [8]).This method gives a relationship between the level of noise in the traveltime measures and the estimation error of the velocity distribution.Such alternative enables us to find the parameters of the velocity that are identifiable by inversion together with the traveltimes measures which are significantly involved in the inversion.This sensitivity analysis method is compared with the classical resolution matrix method which is usually used in the seismic inversion.To present the advantages of SVD based method over resolution matrix one, we propose a study case taken from a seismic survey.The choice of our study case is motivated by two points.First, the classical resolution matrix method is not appliable in this case but the SVD based one gives reliable results.Second, we are considering a simple and realistic case of seismic inversion where it is easy to predict intuitively the sensitivity analysis results and to verify the validity of the method.

THE FORWARD PROBLEM
The forward problem seeks to compute the traveltime of a seismic wave knowing the velocity distribution.The PDE that relates these two quantities is the Eikonal equation.This equation represents the basis of our work.It is shown here how it can be derived from the wave equation.For simplicity, the wave equation is considered only with acoustic case (seismic waves in a fluid).The derivation of the Eikonal equation from elastic wave equation (seismic waves in a solid) is very similar but it takes more steps in the mathematical development.The geological domain Ω that represents the underground where the seismic survey is done is two-dimensional and consists of two unbounded horizontally supperposed layers (Ω 1 , Ω 2 ) separated by an unknown interface γ (see Figure 1).

THE EIKONAL EQUATION
Eikonal equation is derived from the waves equation: (2) where π : Ω ϫ R + → R is the pressure field due to the wave's propagation and υ(x) is the spatially varying seismic wave's velocity.Practically, such wave will be referred to as "seismic ray".Its main features are as follows (see 1): 1. is defined on a compact support subset of Ω that includes the source, the receiver and also intersect γ.The velocity of the seismic wave in the domain is supposed to be continuous except along the interface γ. 2. Thus, ray shooting must have an incidence angle that allows fulfillment of point one.
The Eikonal approach to the Eqn.(2) seeks solutions in the high frequency domain, where an appropriate approximation of the solution π is sought in the form: (3) The functions Π(x) and T (x) represent respectively the amplitude and the traveltime of the wave at point x.
Injecting equation ( 3) in (2) leads to: By neglecting ٌ 2 Π = 0 and by setting all the coefficients of the powers of ω to zero, one deduces two main equations, the first of which is the transport equation obtained by setting to zero the coefficient of ω : ( The second equation obtained by setting to zero the term of Eqn.(4) in ω 2 gives Eikonal model: (6) This is a nonlinear first order PDE relating the wave's traveltime T(x) to the wave's velocity.If we let be the components of the slowness vector and be the Hamiltonian function an equivalent form of the Eikonal equation is: More details details about the Eikonal equation or the seismic ray theory in general can be found in [9] where the existence of the solution is obtained using the method of characteristics.This method seeks to replace the nonlinear first order partial differential equation (7) with a system of ordinary differential equations along a special curve called characteristic or seismic ray travel curve.Specifically, let x(s) be a parametric representation of the characteristic curve with a parameter s, T (s) = T (x(s)) and p(s) = ٌT (x(s)).By by formally differentiating equation (7) respect to s, we obtain: This equation is satisfied when the characteristics obey to the follows system (see( [10])): (9) Differentiating T (x(s)) along the characteristic curve we obtain an additional information for transporting T along the curve: (10) Replacing H with 1  2 {||p|| 2 -1/υ 2 }, Eikonal equation can be equivalently formulated through the characteristic system: (11) where s ʦ [0, s max ] with s max to be determined by solving the first four equations of (11) (that are independent from T ) in a way whereby the seismic ray originates at the source, intersects (γ) and ends up by reaching the receiver when s = s max .The resulting travel curve (Γ), is then used to obtain the the total traveltime by integrating the fifth equation over (Γ): (12) The initial conditions on x(s) and p(s) are given by (13) ( ) where α is the direction of the ray at the source, to be determined by the shooting method technique.

THE SHOOTING METHOD
In the case of the seismic survey, each ray is composed of three segments (see Figure 1) and has a discontinuity when the ray hits the interface γ.The discontinuity is handled using the Snell-Descartes law (see [9]).The computation of the ray is done in three steps: 1.The ray leaves the source (x(0) = x source ) with an angle α with respect to the horizontal such that α > α 0 ; so this minimal angle α 0 would permit the ray to hit the interface γ and penetrate in the second layer.As such, the ray direction is given by: p(0) = (cos(α)/υ(x(0)), sin(α)/υ(x(0))).
The initial value problem ( 11) is then solved in Ω 1 using improved Euler scheme leading to an intersection point Q 1 = x(s 1 ), s 1 being the curvilinear coordinate of Q 1 .2. The Snell-Descartes law is then used to obtain p 1 , the new direction of the ray in Ω 2 , from (x(s 1 ), p(s 1 )) in Ω 1 .Similarly to step 1, the initial value problem (11) is solved in Ω 2 until the ray hits the interface at Q 2 = x(s 2 ), the second point of intersection of the ray with the interface which curvilinear coordinate s 2 .3. Once more the new direction of the ray in Ω 1 , p 2 , is obtained from (x(s 2 ), p(s 2 )) using Snell-Descartes law and System (11) is solved in Ω 1 with in initial conditions (x(s 2 ) = Q 2 and p(s 2 ) = p 2 ) until the ray hits the surface of the domain at point x final .Note that x final depends on α, the direction of the ray at the source.The shooting method technique aims to find α for which x final (α) = x receiver so that the computed ray connects the source to the receiver.The implementation of the method considered is in two steps: 1.
Define for α an interval of existence [α min , α max ] such that both α min and α max > α 0 .In our study case, one must have |α min -α max | ∼ ∼ 10 −3 rad and x receiver belongs to the segment defined by x final (α min ) and x final (α max ).

2.
Shoot N rays having α ʦ [α min , α max ], N ∼ ∼ 1000, then pick the two values of α, α 1 and α 2 for which | x final (α)x receiver | is the smallest.Given that the width of the interval [α min , α max ] is O(10 -3 ), this implies that α will be determined with at least 6 significant figures and x f inal (α) is very sensitive to α.
However, for some velocity distributions (υ 1 , υ 2 , γ), solving the equation x final (α) = x receiver may require N shots of rays with N > 1000 and the simulation of the forward problem becomes very expensive in computations.To overcome this difficulty, we propose in section 5 a formulation that avoids such expensive ray shooting thus allowing the inversion to become feasible in time execution.

MODEL SPACE AND DATA SPACE
The model space, the space in which we are seeking the velocity distribution (υ 1 , υ 2 , γ), has to guarantee the existence of G the Jacobian of BT (υ) because it will be needed to perform the p x x ( ) (cos( )/ ( ( )), sin( )/ ( ( ))) 0 0 0 sensitivity analysis.The theoretical basis of the choice of the functional space of (υ 1 , υ 2 , γ) is given in ( [11]).One must seek for (υ 1 , υ 2 , γ) in = (H 1 (Ω 1 ), H 1 (Ω 2 ), H 1 (Γ)) with H 1 the Sobolev first order space .For practical reasons has to be discretized into a finite space using a regular grid of points and some interpolation functions.is the space of the unknown velocity and R R + NT is the data space, where N T is the number of traveltimes measures.
Bilinear function's are used to define .We consider to represent the domain with a rectangular grid of three horizontal lines and M vertical lines (see Fig. 2).The coordinates of the nodes of the grids are in {0, z 1 , z 2 } ϫ {0, x 1 , . . ., x M } where {0, z 1 , z 2 } are the depths of the horizontal lines and {0, x 1 , . . ., x M } are the positions of the vertical lines.We denote by υ 2 is replaced by a bilinear function that interpolates it at the four nodes of the cell.Thus in each cell υ(υ) is approximated by υ h (x) such that: (15) The discretisation of the velocity in the upper layer υ 1 1 h depends on the the values of υ 1 at the nodes of the grid in the upper layer: (υ 1 (x 1 , 0), . .., υ 1 (x M , 0)) and (υ 1 (x 1+M , z 1 ), . .., υ 1 (x 2M , z 1 )).
The interface is taken as a broken line having its nodes on some edges of the grid.The coordinates of the nodes of the interface are (x i , z i int ), i = 1, . .., M where {z i int } are the depths of the nodes of the interface.
We denote by m the vector of all the parameters on which depend the discreet velocity distribution (υ 1 h , υ 2 h , γ h ): ), ( ), . .., ( , ) , Parametrization of a geological domain of two layers: bilinear interpolation of in each cell and a broken line interface.
is defined as the set of functions υ(x) such that 1/υ 2 is a bilinear function over each cell of the grid and discontinuous along the interface between the two layers.The inverse problem is seeking m ʦ h given the array of traveltime measures d R R + NT .

FRÉCHET DERIVATIVE OF THE TRAVELTIME WITH RESPECT TO THE VELOCITY
Our study is restricted to the sensitivity of the traveltimes with respect to the velocity model only since it is known that in the wide angle context traveltimes are less sensitive to interface's shape than velocity (see [4]).On the basis that the pair (T, υ) solves the inverse problem T = F (υ), where F is the eikonal solver, we seek the Fréchet derivative F' (υ) so as For such purpose we start by linearizing Eikonal equation.If Γ k is the ray connecting a fixed source to a receiver k, let: where υ(x) is the velocity field associated with the seismic ray Γ k .Let now δ be the perturbation operator: δυ and δT k are the perturbations of the velocity and the traveltime respectively.A first order Taylor expansion about υ leads to: The integral over δΓ is equal to zero from the Fermat principle with fixed endpoints (see [4]).Thus we find: On the other hand note that Γ k,i being the intersection of Γ k with a grid cell i.Replacing υ by υ h and using expression (15): Assembling (23) for all the cells i of the grid and all rays k, the linearized Eikonal equation becomes in a matrix form: where, a ʦ R 8M is an array of the coefficients of all the bilinear forms used to approximate over the grid cells and J ʦ R N T × 8M the Jacobian matrix deduced from Eqn. (23).
Equation ( 25) can be rewritten in terms of the value of at the grid nodes whereby a = Cu with C ʦ R 4M؋8M being the matrix of bilinear interpolation.Hence: (26) and is the Jacobian matrix we are seeking to carry out the sensitivity analysis.

INVERSE FORMULATION
A classical formulation of the inverse problem would be: In this paper we do not give theoretical results about the existence and the uniqueness of the solution of this optimization problem (see [12] and [13]).However we know from ( [11] and [14]) that this optimization problem is ill posed and one has to add regularization terms and/or constrains to the cost function in order to guarantee the well posedness of the problem.In the following two subsections we present the regularized version and the constrained version.

REGULARIZED INVERSE FORMULATION AND RESOLUTION MATRIX BASED SENSITIVITY ANALYSIS
With the first possibility, the inverse problem is formulated as follows: (28) where m 0 is some reference discreet velocity distribution, C m is a covariance matrix of the discreet velocity distribution and η is a damping factor.The solution of (28) is the discreet velocity distribution m * such that the gradient of the cost function in ( 28) is null at m * : with G is the jacobian matrix This equation can be solved using an iterative method such as the Gauss-Newton method.A sequence of discreet velocity distributions m k is constructed as follows: (30) We keep updating m k until the cost function is less then a tolerance value.
Adding and retrieving to Eqn.This linear relationship between m * and m true stats that each element of m * is a linear combination of the elements of m true .The sensitivity analysis based on the resolution matrix consists in checking the values of the diagonal elements of R. If R ii is close to 1 then the element m * i is estimated with a high precision.If R ii < = 0.7 then it is usually considered that the estimation of m * i is not reliable.The major drawbacks of this method are that: • it requires providing the covariance matrix C m of the model.If this matrix is not provided (C m = 0), this leads to R = I and m * = m true which is not possible using noisy measures.So this method is only appliable in the case of the regularized formulation and not in case of the constrained formulation, • it assumes that ξ ≡ 0 so there is no direct relationship between the measure errors and the accuracy of the estimated velocity distribution.With this method, it is not possible to know which parameters of the discreet velocity are not identifiable for a certain level of noise in the traveltimes measures.

CONSTRAINED INVERSE FORMULATION AND SVD BASED SENSITIVITY ANALYSIS
We choose to add constrains to 27 based on physical considerations.As outlined above, the refracted rays must have a certain curvature in order to be able to connect the source with the receiver, and given that the curvature of the ray is proportional to the gradient of the velocity we impose constrains on the gradient of the velocity to insure the wellposedness of the problem.The gradient of the velocity in the first and the second layer is considered bounded between the values G 1, min , G 1, max , G 2, min , G 2, max respectively.We write also the implicit ( ) constrain that the number of computed traveltimes N R (m) for the discreet model m is equal to the number of traveltimes measures N T . (33) To solve 33, one can use an iterative deterministic method such as the SQP augmented Lagrangian method in [14].This method constructs a sequence of iterations m k until iteration k such that , the mean square misfit error, is less then a tolerance ε.One can note that the choice of ε is independent of N T , the number of traveltimes measures.The difficulty in our case is that the computation of all the elements of the vector BT (m k ) is not always possible.
Let B c (m) be the operator that gives the available computed traveltimes, and A(m)d the corresponding traveltimes measures, where A = A(m) ʦ R N R (m) × N T is the correspondence matrix ((A) ij (m) = 0 or 1) between the measured data and the associated simulated traveltimes.
is the mean square misfit error for the N R (m) computed traveltimes.Our formulation of the inverse problem consists in finding the discreet velocity distribution m for which the number of computed traveltimes is maximal under the constrain that (34) We consider that formulations (1) and ( 2) are equivalent in the sens that: • it is obvious that if m * 1 is a solution for formulation 1 it is also a solution for the formulation 2: m * 1 is the solution of formulation 2 which gives 2 is a solution for the formulation 2 then We found that the constrained version gives more reliable results in our case, so we will be using this method in the numerical experiment.
We carry out a sensitivity analysis following the one given in [8].
Let n = N T , q = 4M (n >> q) and consider the singular value decomposition (SVD) of the Jacobian matrix (F'(u)) which was defined in section (4).
where W = [w 1 ,...w n ] and V = [υ 1 , ..., υ q ] are orthogonal basis for traveltimes and velocity respectively: One easily verifies that: Given that A first order analysis yields to: Hence we have: Otherwise, if there exists some k such that then the estimation of the k th parameter in the singular basis can't be validated ([8]).Inequality (38) represents the forward relationship between the traveltime measure error and the error in the velocity estimation.
It will be used in the seismic inversion in section 7.

DATA ACQUISITION
In 2002, the IFREMER 1 has done a seismic survey in Morocco.The goal of this survey is to study the structure of the margin of the Moroccan coast in the sea and in the underground.Figure (3) shows the geographical location of the survey.A linear profile 50km long is considered, 28 shoots (red stars) were done at the surface of the ground by shacking trucks.The signals were recorded by 25 different stations (black dots).
The recorded signals are processed to eliminate the noise as well as any other perturbing factor.The signal process is in four steps: 1. Correlation: the recorded signal is compared with the source signal to extract the significant response of the underground.In our case study, the observation of the seismic sections shows the existence of two signals which suggests the existence of two supperposed layers separated by an interface.This is the reason of considering in this paper a two dimensional domain of two layers.Typically, these two layers would be a layer of basalt (υ 2 Ϸ 5000m.s−1 ) over a layer of sediments (υ 1 Ϸ 2000m.s - ) .One can deduce from these sections the traveltimes of the first arriving signals.Traveltimes represents the input data of the inverse problem.They are of type wide angle refraction because the first arriving signals have to propagate in the second layer where the velocity is much higher than in the first layer.This hypothesis with the range of distances between the source and the receiver put some constrains on the shape of the refrated rays.The ray must have an incident angle with the interface very close to the critical angle and it has to be slightly bended in the second layer to be able to propagate over a long distance near the subsurface (see Figure 1).

NUMERICAL RESULTS
Traveltime measures are taken from the seismic experiment described in section 6 with N T = 102 and ε = 0.1s.G 1, min , G 1, max , G 2, min and G 2, max are defined on the basis of prior geophysical information on the nature of underground.We solve the inverse problem with different values of M, the number of columns in the grid, to determine the number of columns needed to model accurately the physical domain.This approach is called a multi-scale approach.Results for different values of M are shown in table (1).One can deduce that the solution model with M = 5 (N R (m * ) = 76) is sufficient to describe the underground structure studied in the survey.Velocity distribution and ray tracing corresponding to the solution model m * are shown in Figures (5) and (6).
Next, a SVD of the Jacobian matrix about the solution m * with M = 5 is computed.Figure (7) shows the singular values ratios (s k /s 1 , k = 1 . . .q).Following the criteria given in Eqn.(38) to define an identifiable parameter, one can deduce from Figure (7) that four singular discreet velocity parameters (discreet velocity parameters expressed in the singular basis) can be recovered if the level of noise in the traveltimes measures is less than 10%.Given that the absolute traveltime measure error ε = 0.1, we found that among the N T traveltime measures available from the experiment 27 present more than 10% of noise.Thus these measures must be discarded from the inversion to identify correctly four singular parameters.The inversion is repeated with this new subset of traveltime measures (N T = 75) and the discreet velocity distribution obtained now is the "final" discreet velocity distribution containing a reliable estimation of the four singular parameters.Figure (8) shows the "'final" discreet velocity distribution corresponding to the new traveltime measures subset of the survey and Figure (9) shows the corresponding ray tracing.
Figure (10) shows the coordinates of the four singular parameters of the discreet velocity distribution, parameters expressed in the singular basis.Each singular parameters is a linear combination of the physical parameters.The coordinates (1-5, 16-20) of these singular parameters are neglectable when compared to the others.Thus these four singular parameters depend only on the physical parameters numbered (6)(7)(8)(9)(10)(11)(12)(13)(14)(15).Hence it can be deduced that the physical parameters (6)(7)(8)(9)(10)(11)(12)(13)(14)(15) are the identifiable ones.In our numbering, they correspond to the lower nodes of the first layer and the upper nodes of the second layer.Thus one can conclude that with these refraction traveltimes measures, we are able to recover only the velocity near the interface.All the other parameters can not be recovered with the present traveltime measures.
A physical justification would be that, in the considered experiment, traveltimes depend mostly on the second segment of the seismic rays (see Figure (1)).The second ray segment is defined using the velocity in the second layer (the characteristic equations) and velocity in the first layer near the interface (the Snell-Descartes law).
With this technique of sensitivity analysis it is also possible to find the traveltime measures that contain the significant information on the wave's velocity in the domain.From Eqn. (5.2) and given that four singular parameters are identifiable, on can deduce that four singular traveltime (expressed in the singular basis W) constrain effectively the velocity     distribution m.Each singular traveltime is a linear combination of all the physical traveltimes measures.Figure (11) shows the coordiantes of these four singular traveltimes.One can notice that no coordinate is null in all of the four singular traveltimes thus it can be deduced that all the traveltime measures are participating equally in the inversion.

CONCLUSION
In this paper, we presented a traveltime inversion seeking to determine the seismic velocity of shallow layer.The inversion is based on the simulation of seismic rays.In this special case of shallow layer, the computation of some seismic rays may become very in time.
To avoid the computation of such seismic rays, we proposed a formulation of the inverse problem that does not require such expensive rays.
Also we studied the sensitivity of the velocity distribution with respect to the noise in the traveltime measures.The sensitivity analysis is based on a SVD of the jacobian of the traveltimes with respect to the velocity.By relating traveltimes measure errors to the error in the estimation of the velocity, we concluded that we are only able to identify velocity near the interface.A physical justification is given to validate the result.

Figure 1
Figure 1 Seismic ray in a geological domain with two layers.
(29) leads to: (31) Let m true be the true velocity distribution, we have d = BT (m true ) + ξ where ξ is the error measures.If we replace in the RHS of Eqn.(31) m * by m true and assume that ξ ≡ 0 we obtain: (5.2) by In eqn.(36) we obtain:(37) Let the accuracies on u and T be respectively σ u and A singular parameter σ

2 .
Summation: signals coming form the same source-receiver pair are summed to obtain an averaged signal by shoot.3. Filtration: the correlation and averaged signals are filtered again to eliminate the residual noise.

Figure ( 4 )
Figure(4) shows a seismic section which is a compilation of all the processed signals.In our case study, the observation of the seismic sections shows the existence of two signals which suggests the existence of two supperposed layers separated by an interface.This is the reason of considering in this paper a two dimensional domain of two layers.Typically, these two layers would be a layer of basalt (υ 2 Ϸ 5000m.s−1 ) over a layer of

Figure 3
Figure 3 Geographical location of the survey: seismic sources (red stars) and recording stations (black dots).

Figure 4 Figure 5
Figure 4 Seismic section: traveltimes of the first arriving signals.

Figure 6
Figure 6 Seismic traveltime inversion using 102 traveltime measures: the corresponding ray trace.

Figure 9
Figure 9 Seismic traveltime inversion using N T = 75 traveltime measures entailing less than 10% of noise: the corresponding ray trace.

Figure 8 Figure 11
Figure 8 Seismic traveltime inversion: velocity distribution using N T = 75 traveltime measures entailing less than 10% of noise.

Figure 10
Figure 10 Coordinates of first four singular vectors showing the identifiable physical parameters.

Table 1
Solutions of the inverse problem for different values of M used to apply the multiscale approach