Hafedh Ben

Hadj Ali

Stéphane Operto
email: operto@geoazur.obs-vlfr.fr.

Jean Virieux
email: jean.virieux@obs.ujf-grenoble.fr.

Velocity model building by 3D frequency-domain full-waveform inversion of wide-aperture seismic data

We assess 3D frequency-domain acoustic full-waveform inversion data as a tool to develop high-resolution velocity models from low-frequency global-offset data. The inverse problem is based on a classic gradient method. Inversion is applied to few discrete frequencies allowing management of a limited subset of the 3D data volume. The forward problem is solved with a finite-difference frequency-domain method based on a massively parallel direct solver allowing efficient multiple-shot simulations involving several thousands of sources.

The inversion code is fully parallelized for distributed-memory platforms taking advantage of a domain decomposition of the modeled wavefields performed by the direct solver. After validation on simple synthetic tests, full-waveform inversion was applied to two targets 1 (channel and thrust system) of the 3D SEG/EAGE overthrust model corresponding to 3D domains of 7 × 8.8 × 3.3 km 3 and 13.5 × 13.5 × 4.65 km 3 respectively. The maximum inverted frequencies were 15 Hz and 5 Hz for these 2 applications respectively. A maximum of 20 dual core biprocessor nodes with 8 gigabytes of share memory per node was used for the second case study. The main structures were successfully imaged at a resolution scale consistent with the inverted frequencies. This study confirms the feasibility of 3D frequencydomain full-waveform inversion of global-offset data on large distributed-memory platforms to develop high-resolution velocity models. Further work is required to [i] perform more representative applications on larger computational platforms, [ii] assess the sensitivity of the 3D full-waveform inversion to the acquisition geometry and to the starting model and [iii] assess whether velocity models developed by full-waveform inversion can be used as improved background model for prestack depth migration. Phase (radians) Receivers along y Receivers along y Receivers along x Receivers along x Receivers along x Receivers along x Receivers along x Receivers along y Receivers along x Receivers along y Receivers along x Receivers along x Receivers along x Receivers along x

INTRODUCTION

Three-dimensional quantitative seismic imaging in complex environments such those involving deepwater, thrust belts, subsalt and subbasalt structures is one of the main challenge of seismic exploration for hydrocarbon exploitation. In the depth domain, the imaging flowchart for multichannel seismic reflection data is subdivided into 2 main steps: estimation of a velocity macromodel and prestack depth migration (PSDM). The first task is critical since the velocity macromodel has a strong impact on the accuracy of the migrated images in terms of focusing and positioning in depth of the reflectors. The criteria that must verify the velocity macromodel to provide accurate migrated images are still unclear (see [START_REF] Operto | Can we image quantitatively complex models with rays?[END_REF][START_REF] Operto | 3-D ray-Born migration/inversion. part 2: application to the SEG/EAGE overthrust experiment[END_REF]) for some illustrations of the sensitivity of 2D and 3D true-amplitude PSDM to the accuracy of the velocity macromodel). Estimation of a reliable velocity macromodel for PSDM from conventional multichannel seismic reflection data is a difficult task which is even more dramatic in complex environments because of the velocity-depth ambiguity at significant depths. The most commonly-used approaches to build velocity model for PSDM rely either on reflection traveltime tomography (e.g., Stork (1992)) and migration velocity analysis (e.g., [START_REF] Yilmaz | Migration velocity analysis by wavefield extrapolation[END_REF]; [START_REF] Chauris | Migration velocity analysis from locally coherent events in 2d laterally heterogeneous media, part 1[END_REF]). Both approaches make approximations for modeling wave propagation such that the high-frequency approximation or the one-way approximation of the wave equation. This incomplete modeling of wave propagation together with the limitations imposed by narrow-aperture acquisition geometries, can preclude imaging of steeply dipping reflectors.

For example, it was recently illustrated that turning waves and multi-reflected arrivals can contribute to improve imaging of the flank of salt bodies [START_REF] Zhang | Imaging complex salt bodies with turning-wave oneway wave equation: 76th annual meeting[END_REF]. As a result, many efforts has been dedicated these last years to the extension of one-way propagator to exploit these arrivals in PSDM (e.g., [START_REF] Zhang | Explicit marching method for reversetime migration: 77th annual meeting[END_REF]).

In this paper, we investigate frequency-domain full-waveform inversion (FWI) of wideaperture data as a tool to build 3D high-resolution velocity model in complex environments (Pratt, 2004). By wide-aperture acquisition survey, referred in the following as global-offset acquisition, we mean any acquisition system with sufficiently long offset coverage to allow recording of turning waves whose refraction depth cover the zone of interest. Wide-aperture arrivals are primarily sensitive to the large and intermediate wavelengths of the medium which are difficult to image from multichannel seismic reflection acquisition and limitedbandwidth sources. Moreover, multifold wide-aperture surveys lead to a redundant control of frequency and aperture angle on the wavenumber illumination in the model space. This redundancy can be decimated to design efficient numerical approaches for seismic imaging in the frequency domain [START_REF] Pratt | Seismic waveform inversion in the frequency domain, part I : theory and verification in a physic scale model[END_REF][START_REF] Pratt | Inverse theory applied to multi-source cross-hole tomography. part 1 : acoustic wave-equation method[END_REF]. Global-offset acquisition survey may be carried out at sea and on land with a newtwork of ocean bottom and land stations respectively (see [START_REF] Clarke | Processing of a novel deepwater, wide-azimuth node seismic survey: The leading edge[END_REF] for a recent 3D wide-azimuth node survey). FWI refers to imaging method based on the complete solution of the full (two-way) wave equation for the forward problem and on inverse problem theory for the imaging problem (Tarantola, 1987). A model is built by minimization of the misfit between the recorded data and that computed in a starting model. The underlying imaging principle shares some similarities with generalized diffraction tomography (Wu and Töksoz, 1987;[START_REF] Pratt | Gauss-newton and full newton methods in frequency-space seismic waveform inversion[END_REF]: the misfit wavefield is processed as the wavefield scattered by the missing heterogeneities in the starting model. These heterogeneities are processed as a series of closely-spaced diffractors. In virtue of Huygens' principle, image of the perturbation model is built by summation of the elementary images of each diffractors. The main drawbacks of FWI methods are two folds: first, they are very computationally expensive due to the complete resolution of the wave equation for a large number of sources. Second, they lack robustness due to the complexity of the global-offset wavefields and their sensitivity to noise and to the inaccuracies of the starting model. During last decade, it was shown in the 2D case that the frequency-domain formulation of FWI applied to global-offset acquisition was providing a promising framework to mitigate these 2 difficulties. Aim of this paper is to start investigating the extension to 3D of this approach.

The frequency-domain formulation of FWI was originally developed for 2D cross-hole acquisition surveys which involve wide-aperture propagations [START_REF] Pratt | Seismic waveform inversion in the frequency domain, part I : theory and verification in a physic scale model[END_REF][START_REF] Pratt | Inverse theory applied to multi-source cross-hole tomography. part 1 : acoustic wave-equation method[END_REF]. Only few discrete frequencies are required to develop a reliable image of the medium thanks to the wavenumber redundancy provided by multifold wide-aperture geometries. Some guidelines to define the optimal frequency interval for FWI are given in Sirgue and Pratt (2004). This frequency-domain decimation leads to very compact volume of data to be managed which may be critical for 3D applications. Second, the frequency-domain formulation of FWI provides the most natural framework to design a hierarchical multiresolution imaging strategy which helps to manage the non linearity of the inverse problem by proceeding from the low frequencies to the higher ones [START_REF] Pratt | Seismic waveform inversion in the frequency domain, part I : theory and verification in a physic scale model[END_REF][START_REF] Pratt | Gauss-newton and full newton methods in frequency-space seismic waveform inversion[END_REF][START_REF] Pratt | Inverse theory applied to multi-source cross-hole tomography. part 1 : acoustic wave-equation method[END_REF]. These two aspects (inversion of a limited subset of frequencies and successive inversion of increasing frequencies) allow to mitigate the two abovementioned drawbacks of FWI. Application of FWI to 2D real data case studies has been limited to frequencies smaller than 20 Hz [START_REF] Ravaut | Multi-scale imaging of complex structures from multi-fold wide-aperture seismic data by frequency-domain full-wavefield inversions: application to a thrust belt[END_REF][START_REF] Operto | Crustal imaging from multifold ocean bottom seismometers data by frequency-domain full-waveform tomography: application to the eastern nankai trough[END_REF]Shipp and Singh, 2002;[START_REF] Hicks | Reflection waveform inversion using local descent methods: estimating attenuation and velocity over a gas-sand deposit[END_REF]. In 3D, the computational cost of the forward problem suggests that it is difficult to handle frequencies greater than 10 hz for representative problems [START_REF] Operto | 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study[END_REF]. At this resolution scale, the resulting velocity models can be conceived as improved velocity macromodel for PSDM. However, the relevance of FWI velocity models as reference model for PSDM still requires further demonstration. An illustration with a real data case studies is provided in [START_REF] Operto | Quantitative imaging of complex structures from multi-fold wide aperture seismic data[END_REF][START_REF] Operto | Multiple arrival migration/inversion of real global aperture seismic[END_REF] in the frame of the imaging of a thrust belt in the southern Apennines by combined first-arrival traveltime tomography, FWI and true-amplitude ray-based PSDM.

Full wave propagation modeling is a critical issue in FWI methods since it is the most computationally expensive task in the processing flow. In the frequency domain, the forward problem reduces to the resolution of a large sparse system of linear equations for each frequency to be considered. In the 2D case, the few frequencies involved in the inversion can be efficiently modeled for multiple shots using a direct solver [START_REF] Marfurt | Accuracy of finite-difference and finite-elements modeling of the scalar and elastic wave equation[END_REF]. Since the original work of [START_REF] Pratt | Inverse theory applied to multi-source cross-hole tomography. part 1 : acoustic wave-equation method[END_REF], optimal finite-difference stencils were designed for the frequency-domain method based on direct solvers [START_REF] Jo | An optimal 9-point, finite-difference, frequencyspace 2D scalar extrapolator[END_REF][START_REF] Hustedt | Mixed-grid and staggered-grid finite difference methods for frequency domain acoustic wave modelling[END_REF]Stekl and Pratt, 1998). The extension to 3D of this modeling approach was addressed in [START_REF] Operto | 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study[END_REF] who showed that problems of representative size (e.g., the 3D SEG/EAGE overthrust model) can be addressed at low frequencies (< 10 Hz) on currentlyavailable distributed-memory platforms. However, it remains unclear which approach (time versus frequency domain, frequency-domain approach based on direct or iterative solvers) is the most efficient one for 3D full-waveform inversion [START_REF] Plessix | A helmholtz iterative solver for 3d seismic-imaging problems[END_REF]. This may depend on several parameters related to the experimental setup such that the dimensions of the model, the frequency bandwidth, the number of traces in the acquisition and the acquisition geometry. Three-dimensional frequency-domain modeling methods based on iterative solvers have been recently presented in [START_REF] Plessix | A helmholtz iterative solver for 3d seismic-imaging problems[END_REF]; [START_REF] Riyanti | A parallel multigrid-based preconditioner for the 3d heterogeneous high-frequency helmholtz equation[END_REF]Warner et al. (2007). These approaches are far less memory demanding than that based on direct solvers but their run times depend linearly on the number of sources which may be a significant drawback in the case of 3D surveys involving several thousands of shots or receivers.

Aim of this paper is to provide some preliminary insights on the feasibility and rele-vance of 3D frequency-domain FWI for building high-resolution velocity models of isotropic acoustic media. Few applications of 3D frequency-domain FWI to synthetic models have been recently presented by Stekl et al. (2007);[START_REF] Clarke | Processing of a novel deepwater, wide-azimuth node seismic survey: The leading edge[END_REF][START_REF] Ben-Hadj-Ali | 3D acoustic frequencydomain full-waveform inversion[END_REF]. Some critical issues related to 3D FWI will not be addressed in this paper and will be left for further investigations. These issues are the building of a reliable starting model for FWI which is conventionally carried out by first-arrival traveltime tomography when global offset acquisitions are considered (Brenders and Pratt, 2006a,b;[START_REF] Ravaut | Multi-scale imaging of complex structures from multi-fold wide-aperture seismic data by frequency-domain full-wavefield inversions: application to a thrust belt[END_REF]).

We will assume in the following that a starting model describing the long wavelengths of the true medium is available. The second issue is the design of 3D global offset acquisition surveys suitable for 3D first-arrival traveltime tomography and FWI. Some footprints of the azimuth coverage on 3D FWI was recently illustrated by [START_REF] Clarke | Processing of a novel deepwater, wide-azimuth node seismic survey: The leading edge[END_REF] thanks to an application to the 3D SEG/EAGE overthrust model. Numerical examples presented in this paper are focused on surface global-offset surveys carried out with networks of source and receiver on the surface.

In the first part of the paper, we briefly review the theory of frequency-domain fullwaveform modeling and inversion. In the second part, we discuss the parallel implementation of frequency-domain FWI for distributed-memory platforms. In the third part, we present several numerical examples of increasing complexities whose aim is to validate the algorithm, to illustrate the sensitivity of the imaging resolution to the acquisition geometry and to provide some insights on the computational complexity of the approach for representative case studies.

THEORY

Theory of frequency-domain full-waveform modeling and inversion will not be rehashed in this paper since it is now a well-established method for imaging 2D media. The extension to the 3D case closely follow the strategies that has been developed in 2D. Therefore, only a brief review of frequency-domain modeling and inversion will be given here. The reader is referred to [START_REF] Operto | 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study[END_REF] for the method used in this paper for frequency-domain wave propagation modeling, [START_REF] Pratt | Gauss-newton and full newton methods in frequency-space seismic waveform inversion[END_REF] for theoretical aspects of frequency-domain FWI and [START_REF] Pratt | Seismic waveform inversion in the frequency domain, part I : theory and verification in a physic scale model[END_REF] for more practical aspects such that waveform inversion data preprocessing and source estimation.

3D acoustic finite-difference frequency-domain modeling

The 3D visco-acoustic wave equation in the frequency domain is given by ω 2 κ(x,y,z) P (x, y, z, ω) + ∂ ∂x 1 ρ(x,y,z) ∂P (x,y,z,ω) ∂x

+ ∂ ∂y 1 ρ(x,y,z) ∂P (x,y,z,ω) ∂y + ∂ ∂z 1 ρ(x,y,z) ∂P (x,y,z,ω) ∂z = -S (x, y, z, ω) (1)
where ρ (x, y, z) is density, κ (x, y, z) is the bulk modulus, ω is frequency, P (x, y, z, ω) is the pressure field and S (x, y, z, ω) is the source. Sponge-like PML (Perfectly Matched Layer) at edges of the necessary limited numerical model can be easily implemented in the frequency domain to absorb energy [START_REF] Berenger | A perfectly matched layer for absorption of electromagnetic waves[END_REF][START_REF] Operto | 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study[END_REF].

Since the relationship between the pressure wavefield and the source is linear, the discrete acoustic wave equation (1) can be recast in a matrix form as

A p = s (2)
where the complex-value impedance matrix A depends on the frequency and the medium properties.

The system 1 can be discretized with the so-called parsimonious mixed-grid finitedifference method [START_REF] Jo | An optimal 9-point, finite-difference, frequencyspace 2D scalar extrapolator[END_REF][START_REF] Hustedt | Mixed-grid and staggered-grid finite difference methods for frequency domain acoustic wave modelling[END_REF][START_REF] Operto | 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study[END_REF]. The mixed-grid discretization which consists of the discretization of the differential operators on different rotated coordinate systems is complemented by a mass-term distribution (an anti-lumped mass) which allows to significantly improve the accuracy of the stencil [START_REF] Marfurt | Accuracy of finite-difference and finite-elements modeling of the scalar and elastic wave equation[END_REF]. The combined use of the mixed coordinate systems and the mass distribution allows to design both accurate and spatially-compact stencil. A dispersion analysis demonstrates that only 4 grid points per wavelength are enough to obtain accurate simulations in homogeneous media. This discretization rule is optimal for FWI whose resolution limit is λ/2 where λ is the wavelength. Use of compact stencil is critical if a direct method is used to solve the system resulting from the discretization of the Helmholtz equation because compact stencils allow to limit the numerical bandwidth of the matrix and hence its fill-in during LU factorization. Use of direct solver is interesting in the case of multiple-shot simulations as those required by tomographic applications since the LU factorization is independent of the right-hand side terms in equation 2. To solve system 2, we use the massively parallel direct solver MUMPS which has been developed for distributed-memory platform [START_REF] Amestoy | Multifrontal massively parallel solver (MUMPS version 4.6): Users'guide[END_REF]. A detailed complexity analysis of this approach is provided in [START_REF] Operto | 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study[END_REF].

Frequency-domain full-waveform inversion

The inverse problem is solved by a classic weighted least-squares gradient method (Tarantola, 1987). Newton and quasi-Newton (Gauss-Newton) methods were rejected because of the computational cost of either the Hessian or the approximate Hessian [START_REF] Pratt | Gauss-newton and full newton methods in frequency-space seismic waveform inversion[END_REF].

The weighted least-squares cost function is given by

C(m) = ∆d † W d ∆d (3)
where ∆d is the misfit function (the difference between the observed data and the data computed with model m), the superscript † indicates the adjoint (transpose conjugate) and W d is a weighting operator applied to the data which scales the relative contribution of each component of the vector ∆d in the inversion. Minimization of the cost function leads to the following solution for the model perturbation ∆m after scaling and smoothing of the gradient [START_REF] Pratt | Gauss-newton and full newton methods in frequency-space seismic waveform inversion[END_REF][START_REF] Ravaut | Multi-scale imaging of complex structures from multi-fold wide-aperture seismic data by frequency-domain full-wavefield inversions: application to a thrust belt[END_REF][START_REF] Operto | Crustal imaging from multifold ocean bottom seismometers data by frequency-domain full-waveform tomography: application to the eastern nankai trough[END_REF].

∆m i = -α (diagH a + I) -1 G m Re p t ∂A t ∂m i A -1 W d ∆d * . (4)
where diag H a = diagRe{ J t W d J * } denotes the diagonal elements of the weighted approximate Hessian H a and J denotes the sensitivity matrix.

One element of the sensitivity matrix is given by

J k(m,n),i = p t m ∂A t ∂m i A -1 δ n . (5)
where k(m, n) denotes a source-receiver couple of the acquisition system, m and n denotes respectively a shot and a receiver position. δ n is an impulsional source located at the receiver position n.

The diagonal of the approximate Hessian provides a preconditioner of the gradient which properly scales the perturbation model (Shin et al., 2001). The damping parameter is used to avoid numerical instabilities (i.e. division by zero). The matrix G m is a smoothing regularization operator. It is implemented in the form of a 3D Gaussian spatial filter whose correlation lengths are adapted to the inverted frequency component [START_REF] Ravaut | Multi-scale imaging of complex structures from multi-fold wide-aperture seismic data by frequency-domain full-wavefield inversions: application to a thrust belt[END_REF]).

An amplitude gain with offset can applied to each seismic trace within the operator W d

w d (o sr ) = exp(gLog |o sr |) (6)
where the scalar g controls the amplitude of the gain with respect to the source-receiver offset o sr . In our algorithm, the scaling of the gradient could be estimated once per frequency before the first iteration and kept constant over iterations or re-computed at each iteration.

The term ∂A ∂m i is the radiation pattern of the diffraction by the model parameter m i . In the case of the P-wave velocity, this radiation pattern is that of an explosion. In other words, this matrix reduces to one scalar. The source term can be estimated in the FWI algorithm by solving a linear inverse problem [START_REF] Pratt | Seismic waveform inversion in the frequency domain, part I : theory and verification in a physic scale model[END_REF]. The inversion code can be applied to vertical geophone data or to hydrophone data generated by explosive sources.

Indeed, vertical geophone data can be processed as pressure data thanks to the reciprocity principle [START_REF] Operto | Crustal imaging from multifold ocean bottom seismometers data by frequency-domain full-waveform tomography: application to the eastern nankai trough[END_REF]. The inversion is applied in cascade to several groups of discrete frequencies. All the frequencies of one group are inverted simultaneously. The final model obtained close to inversion of one group of frequencies is used as a starting model for the next group of frequencies. For each frequency group, several iterations can be computed.

PARALLEL NUMERICAL IMPLEMENTATION

We use the massively parallel direct solver MUMPS [START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF][START_REF] Amestoy | Multifrontal massively parallel solver (MUMPS version 4.6): Users'guide[END_REF] based on a multifrontal method [START_REF] Duff | The multifrontal solution of indefinite sparse symmetric linear systems[END_REF] to solve the forward problem (system 2). Before LU decomposition, the matrix coefficients are ordered so that dependencies in the graph of the matrix are minimized. Using nested dissection ordering, the theoretical memory complexity of the factorization for a 3D FD problem is O(n 4 ) and the number of floating-point operations is O(n 6 ) where n is the number of grid points along one dimension of the 3D square FD grid (Ashcraft and Liu, 1998). The source vectors for the resolution phase are provided in sparse format on the host processor. After resolution, the multiple solutions are distributed over processors following a domain decomposition driven by the distribution of the LU factors. This means that each processor stores a spatial sub-domain of all the solutions. We take advantage of this distributed in-core storage of the forward problem solutions to solve in parallel the inverse problem.

The central component of the FWI algorithm is the computation of the gradient of the cost function. This operator is basically computed by a weighted summation of the forward problem solutions (FPS), namely, the incident and the backpropagated residual wavefields computed in the starting model, equation 4. The weights in the summation account for the radiation pattern of the diffraction tomography reconstruction (the operator ∂A t ∂m i in equation 5) and for the data residuals. This weighted summation is computed in parallel in a straightforward way by taking advantage of the distribution of the FPSs: each processor computes the subdomain of the gradient corresponding to the subdomain of the FPSs stored on this processor. At the end of the summation, the distributed gradient is gathered on the master processor with a collective communication. Note that, when only the P-wave velocity parameter is involved in the inversion, the matrix ∂A t ∂m i reduces to a scalar located on the i th diagonal. This implies that the gradient at position of m i only depends on the values of the FPSs at this same position. In that case, the parallel computation of the gradient doesn't require any point-to-point communication leading to a parallelism efficiency close to 1 for the gradient computation. Note also that all the FPSs remain in core in the algorithm.

No disk swapping is used in the current algorithm. If no enough memory is available to store in core all the FPSs in addition to the LU factors, the computation of the multi-RHS resolution with MUMPS, of the gradient and of the Hessian is performed in a sequential loop over partitions of RHS terms. Each partition loads in core the maximum number of solutions fitting the available memory. The efficiency of the parallel inversion algorithm is mainly controlled by that of the LU factorization. We obtain a maximum speed-up of 13

with the MUMPS direct solver on our applications [START_REF] Operto | 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study[END_REF].

The parallel FWI algorithm is summarized in Figure 1.

SYNTHETIC EXAMPLES

In this section, we present several numerical examples of 3D FWI of increasing complexities to validate the algorithm and to give some insights on the computational cost of the approach on realistic cases. All applications presented hereafter were computed on a HP DL 145G2 Beowulf cluster of computer center SIGAMM host by Observatoire de la Cote d'Azur (France). This parallel distributed computer is a 48 dual node cluster machine with 2.4 GHz bi-processors with a 19.2 Gflops peak performance per node. This computer has a distributed memory architecture, where each node has 8 GBytes of RAM. The interconnection network between processors is Infiniband 4X. Data sharing among processors is performed using the message passing library MPIHP. For all the examples presented hereafter, the PML layers spread along 5 grid points on each side and each direction. These PML grid points are not taken into account in the description of the finite-difference grids.

3D full-waveform inversion in 2D configuration

In a first step, we validate the 3D FWI algorithm by comparing the results obtained using a 2D FWI code and the 3D one applied in a 2D configuration. Two-dimensional experiments can be designed considering 2.5-D velocity models (laterally invariant in the y-direction)

and an infinite line source in the y direction. The infinite line source in the y direction was implemented on a limited computational domain in the y direction using periodic boundary conditions on the two faces of the model corresponding to y = 0 and y = y max .

The periodic boundary conditions that were implemented are

∂P ∂y y=-h/2,ymax+h/2 = 0 (7)
They are applied on two virtual ghost faces located outside the computational domain at positions y = -h/2 and y = y max + h/2 where h stands for the grid interval.

We applied 3D and 2D FWI to a dip section of the overthrust model [START_REF] Aminzadeh | 3-D Salt and Overthrust models[END_REF] (Figure 2

Inclusion models

In this section, we present application of 3D FWI for simple velocity models composed We first consider the case of a velocity model with one inclusion in the homogeneous background. The velocity in the inclusion is 3500 m/s (Figure 6). The inclusion is centred on the 3D grid. The 4 frequencies were inverted successively. Some horizontal and vertical sections of the inclusion are shown in Figure 6. Note the vertically-elongated shape of the inclusion in the vertical cross-section section and the symmetric shape of the inclusion in the horizontal slice. The vertically-elongated shape of the inclusion is due to the fact the top and bottom parts of the inclusion are mainly sampled by downgoing transmitted wave paths which have a limited resolution power while the shape of the inclusion in a horizontal plane is mainly controlled by reflections associated with shots and receivers located near a same face of the 3D model. This relationship between the aperture illumination and the resolution of the imaging is also illustrated on the two graphs extracted from a vertical and horizontal section running trough the inclusion. The vertical graph exhibits a clear deficit of high wavenumber due to transmission-like reconstruction while the horizontal graph exhibits slight deficit of small wavenumbers due to reflection-like reconstruction. The symmetry of the image of the inclusion in the horizontal plane which results from the symmetry of the inclusion with respect to the acquisition geometry is an additional validation of the 3D FWI algorithm.

The second example contains two spherical inclusions (3500 m/s and 4500 m/s) corresponding to a positive and negative perturbations in the homogeneous background (Figure 7). The centre of the two inclusions lies on the same vertical plane in the middle of the grid. Aim of this test is to verify that the 3D inversion handles properly multi-scattering occurring between the two inclusions. The starting model, the frequencies involved in the inversion and the acquisition geometry are the same that for the previous example. For this case study, the 4 frequencies were inverted both successively and simultaneously (Figures 8 and9). In both cases, the inversion successfully imaged the two inclusions.

Inclusion+interface velocity model

A more realistic example consists of a velocity-gradient layer above a homogeneous layer.

An inclusion was incorporated into the velocity-gradient layer (Figure 10(a)). The minimum and maximum velocities are 3.8 and 6.0 km/s respectively. It is discretized on a 100×100×40 grid with a grid spacing h = 62.5 m, which corresponds to a physical domain of 6.25 km x 6.25 km x 2.5 km. The grid spacing h was kept constant over the successive mono-

SEG/EAGE overthrust model

The 3-D SEG/EAGE overthrust model is a constant density acoustic model covering an area of 20 km x 20 km x 4.65 km [START_REF] Aminzadeh | 3-D Salt and Overthrust models[END_REF]. It is discretized with 25 m cubic cells, representing an uniform mesh of 801 × 801 × 187 nodes. The minimum and maximum velocities in the overthrust model are 2.2 and 6.0 km/s respectively (Figure 12).

Overthrust model: channel target

Due to limited available computer resources, our first application was limited to a small target of the overthrust model centred on a channel. The maximum frequency involved in the inversion was 15 Hz. A horizontal slice and a vertical section of the model are shown in Figure 13. The model covers an area of 7 km x 8.8 km x 2.3 km and is discretized with a grid spacing h = 50 m. This leads to a 141 × 176 × 46 grid. The minimum and maximum velocities are 3.3 and 6.0 km/s respectively. Acquisition system consists of 1452 sources and receivers deployed on the surface according to two nested grids of 44x33 nodes.

The distance between either two sources or receivers is 200 m. We sequentially inverted 5 frequencies ranging from 5 to 15 Hz. For each frequency, we computed 7 iterations. The starting velocity model was obtained by smoothing the true model with a wavenumber filter with a cut-off wavenumber of 1./500 m -1 (Figure 14). The final FWI model provides a lowpass version of the true model (Figure 15). In order to assess the accuracy of the FWI, we low-pass filter the true model in the time domain with a cut-off frequency of 15 Hz to roughly mimic the exact velocity model that would have been inferred by FWI. (Figure 16). To perform this application, we used 60 MPI processes distributed over 15 dual core biprocessor nodes. Each MPI process used 1.5 Gbytes of RAM (see Table 1). Seven iterations of the inversion of one frequency took approximately 45 hours. Table 1 gathers information related to running time and memory requirement for LU factorization, multishot resolutions (both tasks being devoted to the forward problem), gradient and diagonal Hessian computation. Running time for the resolution phase is very small (0.9 s per source)

and illustrates the main advantage of frequency-domain modeling methods based on direct solvers for tomographic applications involving few thousands of source. Computation of the gradient is also negligible in the frequency domain (4 s) thanks to the summation without disk swapping over a very compact volume of data limited to few frequency components.

Moreover, increasing the number of cores in the inversion would have led to a significant reduction of the computational time at the partial expense of memory saving due to memory overhead during parallel factorization [START_REF] Operto | 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study[END_REF].

Overthrust model: thrust target

We now consider the imaging of a significant target of the overthrust model which incorporates the main thrusts of the model (Figure 18). The minimum and maximum velocities are respectively 2.2 and 6.0 km/s. The model covers an area of 13.5 km x 13.5 km x 4.65 km.

Acquisition system consists of two coincident 43 × 43 = 1849 network of sources and receivers deployed on the surface. The distance between either two sources or receivers is 300 m. We sequentially inverted 2 frequencies: 3.5 and 5 Hz. We computed 10 iterations per frequency. For this application, we adapted the grid interval to the inverted frequency. Grid intervals were h=150 m and 100 m for the frequencies 3.5 and 5 Hz. (2007).

For this application, we run 32 processes distributed over 8 dual core biprocessor nodes (4 MPI processes/node) for the 150-m grid (frequency 3.5 Hz) and 60 processes distributed over 20 dual core biprocessor nodes (3 MPI processes/node) for the 100-m grid (frequency 5 Hz) respectively. Note that the number of process per dual core biprocessor node was decreased from 4 to 3 as the size of the problem increases in order to increase the amount of share memory assigned to each processor for large problems. This allows optimization of the memory use at the partial expense of the running time since the memory overhead decreases with the number of process.

The 10 iterations took about 24 hours and 72 hours for the 3.5-Hz and the 5-Hz frequencies respectively. More detailed information are gathered in the table 2.

examples of increasing complexity to validate the algorithm and to give some insights of the feasibility of our approach. Some preliminary applications to the overthrust model suggests that frequency-domain FWI can be successfully applied at low frequencies (< 5 Hz) on limited-size PC clusters to develop 3D velocity models with a maximum resolution of the order of half the wavelength (that is, 400 m for a velocity of 4000 m/s). Successfull modeling was already performed in a significant target of the overthrust model at 7 Hz [START_REF] Operto | 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study[END_REF] using 400 Gbytes of core memory suggesting that full-waveform inversion is feasible at this frequency if large computational ressources are made available. At the resolution scale provided by this range of frequencies (5-7 Hz), the velocity models may be used as improved background models for prestack depth migration.

Aim of our future work is to overcome the memory limitations imposed by the use of a direct solver only. A possible way is to evolve towards hybrid direct-iterative solvers implemented in a domain decomposition method for which the direct solver is applied to subdomains of limited dimension while the iterative solver is applied to the grid points located at the boundaries between the subdomains. The relevance of this approach will have to be demonstrated when a large number of sources is considered.

Analysis of the sensitivity of FWI to the acquisition geometry and to the accuracy of the starting model needs further investigations before being able to move to 3D real data applications.

facility and the support of the staff. We are particularly grateful to A. Miniussi for his help during the installation of the software on the cluster. Finally, this work was carried out within the frame of the SEISCOPE consortium (http://geoazur.unice.fr/SEISCOPE) sponsored by BP, CGG-Veritas, EXXON-Mobil, TOTAL and SHELL.

M EM F ACT O ALL (GBytes) 67 M EM F ACT O P ROC (GBytes) 1.5 T IM E F ACT O (s) 510
T IM E SOLV E ALL (s) 1270
T IM E SOLV E SOU RCE (s) 0.9

T IM E GRADIEN T (s) 4

T IM E diag HESSIANa (s) 3093 Figure 24: Imaging of the thrust system in the overthrust model: Same that for Figure 23 but for a source located in the middle of the receiver array.

Ben Hadj Ali & al -

  (a)), discretized on a 801 × 187 grid with a grid spacing h = 25 m. For the 3D application, the dip section of the overthrust model was duplicated 5 times in the y direction leading to a 3D 801 × 5 × 187 finite-difference grid. A 2D wavefield computed in this 2.5D model with the above mentioned boundary conditions is shown in Figure3.The starting model for inversion is obtained by smoothing the true velocity model with a Gaussian function of horizontal and vertical correlation lengths of 500 meters (Figure2(b)).To obtain stable inversion, the true velocity structure was set in the first 100 meters of the starting model. The 2D acquisition system consists of a line of 200 sources and receivers equally-spaced on the surface. The corresponding 2.5D acquisition system consists of duplicating 5 times the source and receiver lines in the y direction. We inverted sequentially 7 frequencies ranging from 5 to 20 Hz. For each frequency, we compute 10 iterations. The final velocity models inferred from 2D and 3D FWI are shown in Figure4. Some vertical graphs extracted from these models are compared in Figures 5. They are very similar hence, providing a first validation of the 3D FWI algorithm. The agreement between the final FWI models and the true model is also quite good. Some high-amplitude perturbations are still underestimated mainly due to an insufficient number of iterations (Figures 5). These results also give some insights on the high spatial resolution which can be achieved in the velocity models at relatively low frequencies (i.e., < 15 Hz) by FWI of global offset data thanks to the continuous sampling of the wavenumber spectrum up to a maximum wavenumber of 2/λ 15Hz m -1 (λ 15Hz denotes wavelengths corresponding to a frequency of 15 Hz).

  of homogeneous background with one or two inclusions. The models are discretized on a small 31 × 31 × 31 grid with 250 m cubic cells. The velocity in the background medium is 4000 m/s. The inverted frequencies are 1.75, 2.35, 3 and 3.75 Hz. Shots and receivers are uniformly distributed on the top and bottom sides of the 3D model respectively.

  frequency inversions and was set according to the maximum inverted frequency. The starting model for inversion is the velocity-gradient layer extended down to the bottom of the model (Figure 10(b)). Acquisition system consists of two nested networks of 17 × 17 = 289 sources and receivers deployed on the surface. The distance between either two sources or receivers is 312.5 m. We sequentially inverted 5 frequencies ranging from 2 and 12 Hz. We computed 10 iterations per frequency. The final FWI velocity model is shown in Figure 11(a). A vertical graph across the inclusion extracted from the final FWI model is shown in Figure 11(b). The bottom layer is well recovered thanks to the large offset coverage allowing to quantitatively image a broad range of the layer wavelengths. The shape and amplitude of the inclusion ara incompletely recovered with respect to the expected resolution of the imaging at 12 Hz although the image of the inclusion remains clear in Figure 11(a). It may be due to the high velocity perturbation associated to the inclusion (1500 m/s with respect to the background) which makes the linearized inversion difficult. Qualitative inspection of the vertical graph also reveals a slight deficit of small (vertical) wavenumbers in the image of the inclusion and of the bottom layer (this is suggested by the negative velocity perturbations with respect to the true model). This deficit is explained by the surface acquisition geometry which illuminates the vertical components of the wavenumber vector with reflections only.

  Qualitative comparison between the final FWI velocity model and the low-pass filtered true model shows a good agreement between the two models. Comparison between a vertical graph extracted from the starting model, the low-pass true velocity model and the final FWI model is shown in Figure 17. The agreement is reasonably good with again a slight deficit of small wavenumbers in the FWI graph due to the surface-to-surface illumination and underestimation of velocities in the deep part of the model likely due to an insufficient number of iterations.

  These dicretizations lead respectively to grids of dimension 90 × 90 × 32 and 135 × 135 × 47. The starting velocity model was obtained by smoothing the true model with a 3D Gaussian function with a correlation length of 500 meters in the 3 directions (Figure19). The FWI velocity models after inversion of the 3.5-Hz and 5-Hz frequencies are shown in figures 20 and 21 respectively. Although the velocity structure is well imaged at a resolution scale consistent with the inverted frequencies, one can note a checkerboard pattern superimposed near the ends of the horizontal slices of the FWI velocity model. This is clearly a footprint of the coarse acquisition geometry. This pattern is only observed near the end of the horizontal slice because of the non uniform azimuth illumination in these parts of the model. The acquisition footprint has no preferential orientation due to the fact that both shots and receivers are uniformally deployed all over the surface with a constant spacing in the dip and cross directions. Another illustration of the footprint of the azimuth illumination and of the acquisition coarseness on 3D frequency-domain FWI was illustrated by Sirgue et al.

Figure 1 :

 1 Figure 1: Sketch of the FWI algorithm. Tasks performed in parallel are written in green.

Figure 2 :

 2 Figure 2: Imaging of a dip section of the overthrust model: a) True velocity model. b)

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Imaging of a dip section of the overthrust model: (a) Final velocity model from

Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :Figure 13 :Figure 14 :Figure 15 :Figure 16 :Figure 18 :Figure 19 :Figure 20 :Figure 21 :

 78910111314151618192021 Figure 7: Imaging of 2 inclusions by 3D FWI: true model.Ben Hadj Ali & al -

Figure 23 :

 23 Figure 23: Imaging of the thrust system in the overthrust model: a) samplitude (left) and

Table 1 :

 1 Computational cost of the imaging of the overthrust model (channel): M EM F ACT O ALL : Total memory allocated during factorization. M EM F ACT O P ROC : Average allocated memory per working processor during factorization. T IM E F ACT O : Elapsed time for factorization. T IM E SOLV E ALL : Total elapsed time for multi-shot resolution. T IM E SOLV E SOU RCE : Elapsed time for resolution for 1 source. T IM E GRADIEN T : Elapsed time for gradient computation. T IM E diag HESSIANa : Elapsed time for diagonal Hessian

	computation.

Table 2 :

 2 Computational cost of the imaging of the overthrust model (thrust system): F REQU EN CY : inverted frequency. GRID: dimension of the 3D FD grid. M EM F ACT O ALL : Total memory allocated during factorization. M EM F ACT O P ROC : Average allocated memory per working processor during factorization. T IM E F ACT O : Elapsed time for factorization. T IM E SOLV E ALL : Total elapsed time for multi-shot resolution. T IM E SOLV E SOU RCE : Elapsed time for resolution of 1 source. T IM E GRADIEN T : Elapsed time for gradient computation. T IM E diag HESSIANa : Elapsed time for diagonal Hessian computation. Shin, C., K. Yoon, K. J. Marfurt, K. Park, D. Yang, H. Y. Lim, S. Chung, and S. Shin,

	2001, Efficient calculation of a partial derivative wavefield using reciprocity for seismic
	imaging and inversion: Geophysics, 66, 1856-1863.
	Shipp, R. M. and S. C. Singh, 2002, Two-dimensional full wavefield inversion of wide-
	aperture marine seismic streamer data: Geophysical Journal International, 151, 325-344.
	Sirgue, L., J. Etgen, and U. Albertin, 2007, 3D full waveform inversion: Wide versus
	narrow azimuth acquisitions: 77 th SEG Annual Meeting (San Antonio, 23-28 September),
	Expanded Abstracts, 1760-1764, Society of Exploration Geophysicists.
	Sirgue, L. and R. G. Pratt, 2004, Efficient waveform inversion and imaging : a strategy for
	selecting temporal frequencies: Geophysics, 69, 231-248.
	Stekl, I. and R. G. Pratt, 1998, Accurate viscoelastic modeling by frequency-domain finite
	difference using rotated operators: Geophysics, 63, 1779-1794.
	Stekl, I., M. R. Warner, and A. P. Umpleby, 2007, 3D frequency domain waveform inversion
	-synthetic shallow channel example: EAGE 69th Conference and exhibition, Expanded
	Abstracts, C026, EAGE.
	Stork, C., 1992, Reflection tomography in the post migrated domain: Geophysics, 57, 680-
	692.
	Tarantola, A., 1987, Inverse problem theory: methods for data fitting and model parameter
	estimation: Elsevier.
	Warner, M., I. Stekl, and A. Umpleby, 2007, Full wavefield seismoc tomography -iterative
	forward modelling in 3d: 69th EAGE Conference and Exhibition incorporating SPE
	EUROPEC 2007 ,London, Expanded Abstracts, C025, EAGE.
	Wu, R.-S. and M. N. Töksoz, 1987, Diffraction tomography and multisource holography
	applied to seismic imaging: Geophysics, 52, 11-25.

ACKNOWLEDGEMENTS

Access to the high performance computing facilities of MESOCENTRE SIGAMM computer center provided the required computer resources and we gratefully acknowledge both this

We considered that we obtained convergence when the cost function reached about 10 % of its initial value (see Figure 22). We notice that we reached a better convergence rate that shown in Figure 7a of [START_REF] Clarke | Processing of a novel deepwater, wide-azimuth node seismic survey: The leading edge[END_REF]. This may be due to different scaling of the gradient or to the smaller size of our application.

The data fit is illustrated in the frequency domain for the 3.5-Hz and 5-Hz frequencies for 2 shots in Figures 23 and24 respectively. We compare the spectral amplitude and the phase of the monochromatic wavefields at the receiver positions computed in the true velocity model and in the FWI models at the first and last iterations of the two monofrequency inversions. One shot is located at the upper-left corner of the receiver plane (Figure 23) while the second shot is on the middle of the receiver array (Figure 24). The misfit reduction between the first and last iterations is obvious. We note also that this misfit reduction is more effective for the shot located in the middle of the receiver array illustrating the difficulty to match the lower-amplitude arrivals recorded at larger offsets.

CONCLUSION

In this paper, we presented a 3D massively parallel frequency-domain full-waveform inversion algorithm based on a direct solver. Advantages of our approach is related to the robustness of the forward problem provided by the use of a high-performance direct solver, its efficiency to perform multiple-shot simulations in 3D finite-difference grids of relatively small dimensions and a straightforward parallelization of the inverse problem resulting from a domain decomposition of the monochromatic wavefields performed by the direct solver.

Its main drawback is the memory and CPU time complexity of the LU factorization phase which limits the size of the models and the frequency bandwidth which can be addressed on realistic distributed-memory platforms. We presented several applications on synthetic

LIST OF FIGURES

1 Sketch of the FWI algorithm. Tasks performed in parallel are written in green.

RHS stands for right-hand side terms and correspond to sources in the frame of fullwaveform modeling. N proc stands for the number of processes involved in the parallel execution. P i stands for the processor i where P 0 is the master processor. Note that an arbitrary number N of frequencies can be inverted simultaneously (set N f reqgroup = 1 and nf req = N ) or successively (set N f reqgroup = N and nf req = 1). nitermax stands for the maximum number of iterations of one frequency group inversion.


Matthieu Delost

Jean Virieux

Stéphane Operto

First-arrival Travel-Time Tomography using Second Generation Wavelets

Wavelet decomposition of the slowness model has been proposed as a multiscale strategy for seismic first-arrival time tomography. We propose the introduction of so-called second generation wavelets which could be used for any mesh structure and does not need a number of samples as the power of two in each direction. Moreover, one can handle easily boundary effects. A linearized procedure for inverting delayed travel-times considering either slowness coefficients or wavelet coefficients. The ray tracing is solved at each iteration through an eikonal solver while the linear system to be solved at each iteration goes through an iterative solver as LSQR algorithm. We develop wavelet decomposition over constant patches (Haar wavelet) or over linear patches (Battle-Lemarie wavelet) of coefficients at different scales. This decomposition is introduced in the linear system to be solved and wavelet coefficients are considered as unknowns to be inverted. Synthetic examples show that the inversion behaves in a better way as wavelet decomposition seems to act as a preconditioner of the linear system. Local discretisation is possible but requires additional implementation as artefacs once built inside the model description never disappear because of the linearized approach. A binary mask operator is designed for each scale grid and could be applied locally leading to quite different spatial resolution depending on the analysis we could perform of the expected resolution at a given position of the medium. We show that indeed it is possible to design this binary operator and we apply it to synthetic examples as a crosswell experiment inside the Marmousi model. An application to a surface-surface experiment has been performed and the waveled decomposition shows that indeed we may recover detailed features nearby the free surface while preventing imprints of ray coverage at greater depths giving us smooth features at that depths. In spite of the increase demand of computer resources, the wavelet decomposition seems to be a rather promising alternative for controlling the resolution variation of seismic first-arrival tomography.

introduction

First-arrival delayed traveltime tomography (FADTT) applied to wide aperture seismic data has proved to be an interesting tool to investigate the Earth structures [START_REF] Aki | Determination of three-dimentional velocity anomalies under a seismic array using first P-arrival times from local earthquakes[END_REF][START_REF] Spakman | Imaging algorithms, accuracy and resolution[END_REF][START_REF] Hole | Nonlinear high-resolution three-dimentional seismic travel time tomography[END_REF][START_REF] Zelt | Seismic traveltime inversion for 2-D crustal velocity structure[END_REF][START_REF] Benz | Three-dimentional P and S wave velocity structure of Redoubt Volcano, Alaska[END_REF][START_REF] Zhang | Nonlinear refraction and reflection travel time tomography[END_REF]. The resolution of FADTT is controlled both by the theoretical resolution imposed by the relation between velocity structure and first-arrival traveltimes and by the data acquisition configuration as well as the structure itself, which may lead to uneven ray coverage. The theoretical resolution power of FADTT is limited by the size of the first Fresnel zone for each ray. There is an antinomy between the limited resolution power of FADTT and the numerical sensitivity kernel of FADTT corresponding to rays. Rays have no thickness and, therefore, have an intrinsic infinite sensitivity although in practise discretization of the model space always smoothes out this feature [START_REF] Dahlen | Resolution limit of traveltime tomography[END_REF]. Indeed, one raypath induces high frequency information along the ray, even at locations where resolution of FADTT is expected to be low. Inversion can build up models which have a high degree of data fitting but which are far from the true model containing artefacts. Accounting for the limited resolution of FADTT goes through regularized tomographic inversion [START_REF] Farra | Non-linear reflection tomography[END_REF], Delprat-Jannaud & Lailly, 1993[START_REF] Zhang | Nonlinear refraction and reflection travel time tomography[END_REF] often based on the Tikhonov regularization approach [START_REF] Tikhonov | Solution of ill-posed problems[END_REF].

The introduction of the covariance matrix for the model as well as for the data [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least square criterion[END_REF][START_REF] Tarantola | Inverse problem theory: methods for data fitting and model parameter estimation[END_REF] has introduced a precise theoretical background although estimation of those matrices is still a difficult task [START_REF] Monteiller | An efficient algorithm for double-difference tomography and location in heterogeneous media, with an application to kilauea volcano[END_REF]. Another alternative regularization is based on smoothing constraints easier to implement as the laplacian [START_REF] Menke | Geophysical Data Analysis: Discrete Inverse Theory[END_REF][START_REF] Zhang | Nonlinear refraction and reflection travel time tomography[END_REF]. Whatever we choose, regularization strategies lack spatial adaptivity. Frequently, when the smoothing is tuned to remove the footprints of raypaths in poorly illuminated areas, short-scale features that could be resolved in other areas of the model according to the resolution of FADTT are smoothed simultaneously, leading to a loss of information in the well-resolved areas.

The resolution of FADTT generally varies locally. For example, considering surface acquisition geometries, the shallow structures will be sampled both by short-and long-offset rays running across the heterogeneities with different azimuths while the deeper structures will be sampled by long-offset rays only with a much narrower range of azimuths. The resolution of FADTT decreases quite significantly with respect to depth for this geometry. Two other factors creating non-uniform ray coverage are the presence of low velocity anomalies which induces shadow zones in the ray coverage and non-uniform source and/or receiver arrays especially in the case of passive tomography.

In order to consider this spatially varying resolution power, one could introduce adaptive parametrizations based on non-structured grids [START_REF] Michelini | An adaptative-grid formalism for traveltime tomography[END_REF], Spakman & Bijwaard, 1998, Spakman et al., 1998[START_REF] Bohm | 3d adaptative tomography using delauney triangles and vorono 'ı polygons[END_REF][START_REF] Trinks | Adaptative traveltime tomography of densely sampled seismic data[END_REF][START_REF] Zhang | Adaptative mesh seismic tomography based on tetrahedral and voronoi diagrams: application to parkfield, california[END_REF]. The size of the elementary cells are adapted to the local ray coverage such that the amount of control for each model parameter tends to be uniform. More promizing is the multigrid approach where the inversion procedure recovers values on different grids depending on the resolution power [START_REF] Madych | Tomography, approximate reconstruction, and continuous wavelet transforms[END_REF][START_REF] Chiao | Multiscale seismic tomography[END_REF][START_REF] Chiao | Multiresolution parametrization for geophysical inverse problems[END_REF][START_REF] Zhou | Multiscale traveltime tomography[END_REF] based on wavelet tools [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. In this article, we shall propose an further investigation of this multigrid approach using new wavelet transform named as second-generation wavelets [START_REF] Sweldens | The lifting scheme, a construction of second generation wavelets[END_REF].

Wavelet transformation

Over the last years, many constructions of wavelets have been introduced both in mathematical analysis and in the signal processing literature. In mathematical analysis, wavelets were originally constructed to analyse and describe geophysical signals among other signals using translations and dilations of one specified function called mother wavelet. A mathematical framework has been set up by the "French school" [START_REF] Daubechie | Painless nonorthogonal expansions[END_REF][START_REF] Grossmann | Decomposition of hardy functions into square integrable wavelets of constant shape[END_REF]. In signal processing, wavelets found their ways in the context of quadrature mirror filters [START_REF] Mintzer | Filter for distorsion-free multirate filter banks[END_REF][START_REF] Vetterli | Filter banks allowing prefect reconstruction[END_REF][START_REF] Nguyen | Two-channnel perfect reconstruction fir qmf structures which yield linear phase analysis ans synthesis filters[END_REF][START_REF] Vaidyannathan | Theory an design of m-channel maximally decimated quadrature mirror filters with arbitrary m, having perfect reconstruction property[END_REF][START_REF] Woods | Subband coding of images[END_REF]. The connection between the two approaches has been performed through the introduction of multiresolution analysis and the fast wavelet transform by Mallat (1989a,b), [START_REF] Meyer | Ondelettes et Opérateurs, I:Ondelettes, II:Opérateurs de Calderon-Zygmund[END_REF]. An important step has been achieved with the construction of orthogonal, compactly supported wavelets by [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF].

Since then, several generalizations to the biorthogonal or semi-orthogonal cases were presented [START_REF] Cohen | Biorthogonal basis of compactly supported wavelets[END_REF][START_REF] Vetterli | Wavelet and filter banks : theory and design[END_REF]. Biorthogonality allows the contruction of symmetric wavelets and thus the use of linear phase filters.

For seismic tomography, wavelets can be used for description of the model at various scales [START_REF] Chiao | Multiscale seismic tomography[END_REF][START_REF] Chiao | Multiresolution parametrization for geophysical inverse problems[END_REF][START_REF] Zhou | Multiscale traveltime tomography[END_REF]. The slowness model is expanded on a wavelet basis and the tomographic system is solved for the wavelet coefficients. The wavelet transform of the model provides a compact multigrid representation of the model thanks to orthogonal or bi-orthogonal basis functions built by translating (localization property) and dilating (scaling property) a mother wavelet and a scaling function. Those functions are basis of functional spaces generaly noted V i for the scaling functions et W i for the wavelet functions. Those spaces verify certain properties and constitute what is called a multiresolution analysis [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. Each level i corresponds to a grid size. Generally, the resolution of each grid decreases by a factor 2 and the mesh spacing is set accordingly when the level i increases. The starting point is the slowness model on the finest grid of the multiresolution representation (generally V 0 ). The wavelet transform applies a cascade of orthogonal restrictions of the slowness model over different approximation grids. The restrictions of the slowness model on a coarse grid is encoded by the scaling coefficients. Those coefficients are the components along the scaling functions (generally noted φ i ), basis functions of space V i . At each step, the increment of information lost during the restriction on the coarser grid is encoded by the wavelet coefficients before proceeding to the next grid. Those coefficients are associated to wavelet functions (generally noted ψ i ), the basis functions of the spaces W i . At the last iteration of the orthogonal transform, we end up with the wavelet coefficients on each grid plus the scaling coefficients on the coarser grid (see figure 1 for a simple example with a 1D signal). The inverse wavelet transform proceeds in the other direction from the coarse grids to the finer ones.

Multidimensional wavelet basis can be built by a tensor product of 1D basis. The decomposition of the 2D slowness model of a wavelet basis can be written as

u(x, z) = J j=1 +∞ n=-∞ I i=1 +∞ m=-∞ c ww j,n,i,m ψ j,n (x)ψ i,m (z) + J j=1 +∞ n=-∞ +∞ m=-∞ c ws j,n,I,m ψ j,n (x)φ I,m (z) + +∞ n=-∞ I i=1 +∞ m=-∞ c sw J,n,i,m φ J,n (x)ψ i,m (z) + +∞ n=-∞ +∞ m=-∞ c ss n,I,m,J φ J,n (x)φ I,m (z), (1)
where c ww , c ws , c sw and c ss are the wavelet and scaling coefficients. The indices i, j and n, m are the scale and localization indexes respectively. The indices I and J denotes the coarsest levels of the multiresolution approximation. A compact representation of the wavelet coefficients of a 2D signal has been proposed by Mallat (Figure 2).

As Fourier driven transformation, wavelet transformations encounter some relevant drawbacks: the number of samples has to be a power of two, samples have to be regularly spaced and one must deal periodicity which complicates boundary conditions when working on edges or surfaces.

For overcoming these difficulties, the concept of biorthogonal wavelets constructed by lifting scheme has been introduced [START_REF] Sweldens | The lifting scheme, a custom design construction of biorthogonal wavelets[END_REF][START_REF] Sweldens | Buidling your own wavelets at home[END_REF][START_REF] Sweldens | The lifting scheme, a construction of second generation wavelets[END_REF], leading to the so-called second generation wavelets. Wavelets are not necessarily translations and dilations of a defined function. Such wavelets can be adapted to curved edges and surfaces delimited quite complex domains. They are concerned by weighted and irregular samples. The lifting scheme leads to a faster, in-place calculation of the wavelet transform.

In term of signal processing, projections on spaces V and W can be seen as respectively high-pass and lowpass filtering operators. A procedure of direct and inverse transforms can then be interpreted as a perfect reconstruction filter bank (h 1 , g 1 , h 2 , g 2 ), where (h 1 , g 1 ) are the filters for the direct transform and (h 2 , g 2 ) are the filters for the inverse transform. In the orthonormal transform, h 1 = h 2 = h and g 1 = g 2 = g and (h, g) are called conjugate mirror filters. In the biorthogonal case, relaxing the orthogonality condition implies working with larger spaces as Ṽ and W spaces, respectively the dual spaces of V and W . Dual spaces are used for the direct transform and normal spaces are used for the inverse transform (or vice-et-versa). The associated perfect reconstruction filter bank is denoted (h, g, h, g). The lifting scheme is based on the very simple property demonstrated by [START_REF] Sweldens | The lifting scheme, a custom design construction of biorthogonal wavelets[END_REF]:

Given a biorthogonal filter, another one can be built as follows

(h, g, h, g) ⇒ (h, g 1 , h1 , g) H1 (w) = H(w) + S(2w) G(w), G 1 (w) = G(w) -S(2w)H(w),
where capital letters stand for the trigonometric polynomials associated with script letters filters and S is a simple trigonometric polynomial.

This polynomial is chosen so that the wavelet and dual scaling functions satisfy desired properties (as the number Ñ of vanishing moments of the dual wavelet for example). Of course, a dual lifting can also be defined in a similar way:

(h, g, h, g) ⇒ (h 1 , g, h, g1 ) H 1 (w) = H(w) + S(2w)G(w), G1 (w) = G(w) -S(2w) H(w).
With this operation, the real scaling function and the dual wavelet functions can then also be tuned, so they have a requested number N of vanishing moments.

These properties assert that, given a biorthogonal filter, one can construct another one satisfying chosen properties. Starting from a simple biorthogonal filter, it is then possible to bootstrap one's way up to a multiresolution analysis with desired properties on primal and dual wavelets.

The so-called "lifting scheme" proceeds with this strategy: one starts with the straightforward biorthogonal filter that does nothing but set aside odd an even indexed coefficients (the so-called Lazy Wavelet). Then one first does a dual lifting operation by predicting odd indexed coefficients using even indexed ones (we use the spatial correlation between neighbours). We thus ensure that polynomials of degree up to N -1 are well predicted. In a second time, we update the even indexed coefficients using new odd coefficients, so that Ñ moments of the signal are conserved. The algorithm can simply be written as

   (odd i+1 , even i+1 ) = Lazy(λ i ), odd i+1 -= P (even i ), even i+1 + = U (odd i ).
The inverse transform can be written similarly as

   even i+1 -= U (odd i ), odd i+1 + = P (even i ), λ i = M erge(odd i+1 , even i+1 ).
This numerical flowchart can be extended because the predict stage is a simple polynomial interpolation and the update stage is all but the resolution of linear systems: boundary conditions and number of samples are no more a limitation of the wavelet decomposition.

3 Travel-time tomography

Once rays have been computed between each source and each receiver seeing the source through an eikonal solver [START_REF] Podvin | Finite difference computation of traveltimes in very contrasted velocity model : a massively parallel approach and its associated tools[END_REF], the delayed traveltime tomographic system relating slowness anomaly ∆u to traveltime residuals ∆t written as

∆t = A∆u, (2)
can be transformed into the new system ∆t = AW t ∆c,

where the A matrix is the sensitivity or Fréchet matrix and the equation ∆c = W ∆u denotes the transformed slowness model through the transform matrix W . The matrix AW t can be built by computing the 2D wavelet transform of each row of the original matrix A. A row of the sensitivity matrix has n 1 × n 2 number of model parameters in a 2D geometry where the number of horizontal grid points is n 1 and the number of vertical grid points is n 2 : it contains the contribution of one ray to the Fréchet derivatives. After transformation, we have the same number of unknowns but they are in the wavelet domain and partial derivatives are for wavelet coefficients. Once the system equation (3) has been solved through an iterative procedure as the LSQR algorithm [START_REF] Paige | ALGORITHM 583 LSQR : Sparse linear equations and least squares problems[END_REF], the perturbation coefficients ∆u model in the slowness space are obtained by the inverse wavelet transform of ∆c = W ∆u. We shall plot velocity models from these slowness values.

Traveltime tomography in the wavelet domain has an interesting feature: as coefficients are associated to levels of the wavelet transform and, consequently, related to scales and local resolutions, adaptivity strategy is much easier to handle through scale hierarchie. The decomposition of a simple circular anomaly leads to a serie of grids for each scale level using square wavelet basis functions as shown in the figure 3. Controlling incremental perturbations through wavelet coefficients between these different description of the true model taken as the one in the V 0 space is an hierarchical procedure.

As we estimate imprints of rays on different coarser grids through estimations of the Fréchet derivatives, we extend equivalently its area of sensitivity as shown in figure 4. The sensitivity of a ray is related to the spatial discretisation of the slowness model at the given wavelet scale.

The adaptivity of the parameterization can be implemented by applying a mask to each row of the sensitivity matrix AW t , zeroing coefficients in areas of the model where a poor resolution is expected et a given scale.

This adaptivity selection can be implemented via an linear operator M applied to the Fréchet derivatives in the wavelet domain through the following extended linear system ∆t = AM W t ∆c.

(4)

Designing this operator will depend on local resolution analysis. For surface acquisition, one could, for example, define an empirical rule relating the resolution of the model at a given depth to the minimum sourcereceiver offset of the rays sampling this depth. The source-receiver offset would provide an estimate of the width of the first Fresnel zone at a given depth and, hence, of the expected resolution of the tomography varying with depth. According to the sampling theorem, only grids with mesh spacing twice smaller than the width of the first Fresnel zone, would be locally involved in the inversion: a simple rule for construction the M operator. We will discuss later how we can propose a more general strategy for building this operator. Let us remark that the wavelet transform has also the capacity to automatically adapt the parameterization to the resolution in the case where the real resolution is lower than the intrinsic one such as in the case of a shadow zone. In areas where no rays illuminate the medium, partial derivatives of the traveltime time with respect to the wavelet coefficients at the finest scales will be zero, thanks to the second generation property. On the contrary, the partial derivative of the traveltime time with respect to the wavelet coefficients of the coarsest scales will be non-zero due to the spatially extended sensitivity of the coarse grids representations of the rays.

Hence, only the coarse-scale wavelet coefficients of the slowness model will be non-zero after inversion and will automatically lead to a smooth representation of the true model.

Synthetic examples

We present synthetic examples to illustrate how the wavelet decomposition enables us to reconstruct the P-wave velocity model. First, we focus our attention to the wavelet coefficients transformation. We shall show that similar results are found when compared with those of the standard tomography using slowness coefficients.

Finally, we illustrate the efficiency of scales selection using the simple binary operator M . Designing an objective criterion for the spatial definition of this binary operator will be based on resolution matrix analysis.

A densely illuminated model with positive anomalies

As a first example, we design an acquisition system for a complete illumination of the simple circular positive velocity anomaly of 0.5km/s we want to reconstruct (Figure 5) from a background model of constant velocity of 3km/s. This background velocity model will be our starting model.

For each inversion, we estimate the lowest rms error from the misfit function as it flattens with linearized iterations. Moreover, we display horizontal and vertical sections for a better appreciation of the spatial frequency contents.

First, the model is reconstructed using slowness parameters in the inversion over a grid of 5 m step size.

No damping is applied to the over-determined system. The convergence is obtained after four iterations and the anomaly is accurately recovered although with high frequency artefacts (Figure 6). Other artefacts are related with uneven ray sampling in the center of the anomaly. How better will be a reconstruction in the wavelet domain ?

When considering wavelet decomposition, we select wavelets with two vanishing moments. We consider two and four levels for the scale sampling. This means that the coarser grid size will have respectively a 10m resolution and a 40m resolution. The anomaly is well recovered but the convergence to the flat level of the misfit function is slightly slower than before (Figures 7 and8). We need respectively eight and thirteen iterations for reaching convergence. We may underline that the rms final error decreases as we increase the number of levels although the decrease is quite small. For real applications with errors in the picked traveltime, this may not be detected. More importantly, the high frequency artefacts are less noticeable providing a better fit of the model. As we increase the number of levels, we can see the decrease of the amplitude of the high frequency artefacts. In other words, when coarser levels are introduced, the reconstruction is driven by a better spatial coherence, essentially for broad features. We may say that the parametrization into the wavelet domain provides a better conditionning of the system to be inverted, since it is less sensitive to high frequency perturbations.

A densely illuminated model with negative anomalies

We may proceed similarly when considering a negative velocity anomaly. The velocity amplitude will -0.5 km/s and the background velocity model will be the starting model. The reconstruction is more difficult as rays tend to avoid low-velocity zones creating shadow zones which are not sampled by rays (Figure 9).

Velocity anomaly should ne be reconstructed in the center of the medium as rays are missing there. Therefore, reconstructed anomaly comes from the linearized iteration procedure where rays are present in the middle at the early stages of linearized. We may say that values are with constraints and come frome the null space contribution. One must underline that this behaviour is typical of linearized procedure.

The results with the negative velocity anomaly are very similar with those for a positive velocity anomaly, although fondamental differences could be observed (Figure 10). When considering a slowness description of the medium, we observe a poor reconstructed structure in the center. Rays do not sample this zone at the end of the procedure, leaving uncorrected slowness obtained at the initial stages of the inversion. Because of this poor sampling, edge anomalies are observed for correcting this imperfect reconstruction in the middle as few rays of different azimuths are emitted nearby sources and receivers. These edge anomalies are hardly seen for the positive anomaly.

By introducing wavelet coefficients with two vanishing moments, we recover the negative anomaly with an increased number of iterations (Figures 11 and12). The final rms residue is smaller although the decrease is not significant. Because the anomaly we want to recover has a smooth wavelength content, the result is better when considering four levels of discretisation. Still, partial reconstruction introduces artefacts which spread inside the medium from the edges as we move to the four levels description. We can see that wavelet coefficients could not overcome automatically aliasing in the reconstruction essentially based on the linearized procedure. Therefore, we need the operator M for tuning the reconstruction on the specific level we expect in different zones of the medium.

In conclusion of these simple synthetic tests, the wavelet transformation has allowed a better preconditionning of the linearized system to be solved: high frequency perturbations are less noticeable and scale decomposition increases the stability of the reconstructed system at the expense of the number of iterations.

Of course, this analysis is for an ideal coverage by rays and the experiment has been conducted without any noise in the picked times. In spite of this ideal configuration, tuning the parametrization description with the help of expected varying resolution should be added to the automatic wavelet coefficient reconstruction.

A typical surface-surface synthetic example with two different anomalies Moreover, dealing with uneven ray coverage is a crucial issue in tomography: how the wavelet decomposition performs when illumation is far from an ideal one. The experimental data acquisition as well as the structure itself always creates areas of limited ray coverage. Hence parts of the model are expected to be known precisely, while other parts will suffer from low resolution. Wavelets deal with the multiscale nature of the information contained in the data but the inversion is not able to focus automatically anomalies on related scales. Working in the wavelet domain makes easy the selection scales of the parameter model by a simply predefined operator. Zeroing wavelet coefficients corresponding to scales finer than the expected local resolution will help spatial aliasing of the inversion procedure. By applying a mask operator to each row of the sensitivity matrix written in the wavelet domain, we may successfully distinguish well-resolved and poorly-resolved zones in our reconstruction.

Applying a mask does not increase the computational need as we perform a multiplication by zero or one for each term of the matrix. Defining the mask may be quite demanding although pre-tabulated masks should be deviced for typical data acquisition configuration.

The effect of this simple operation is illustrated by the figure 13. In this example, we use two levels of wavelets. the model is a simple constant background model with two heterogeneities a positive and a negative one, enlightened from both sides. For the coarsest level, we keep all the coefficients (the upper left black dial of the mask). For finer levels, we only keep coefficients in a rectangle localized in the center of the model. Every boundary fine coefficients are scaled down to zero. The resulting model is similar to the model constructed whithout mask in the rectangle in the center. On the boundaries, only coarser coefficients have been kept, thus a coarser model coming from the poor ray sampling nearby sources and receivers.

A strongly spatially varying illumination related to a poor ray coverage comes out for a surface acquisition configuration. In this case, the theoretical resolution is expected to be heterogeneous: in the shallow part, both long and short rays spanning across a large range of azimuths make the resolution quite accurate. On the contrary, the deepest parts of the model are sampled and enlightened only by longer rays with a mainly horizontal azimuth leading to low resolution.

Because of the linearized formulation, we thus need a method that can adapt the resolution to the theoretical one aside the natural wavelet reconstruction. Wavelet paramatrization makes it easy to adapt resolution locally, by applying a simple mask. By using an adequate mask, two velocity heterogeneities, a small one in the shallow part and a bigger one in the deep part, are reconstructed according to the expected resolution. The starting model is a constant velocity model without anomalies (Figure 14).

The figure 15 shows the reconstruction through slowness coefficients without any regularization. One can easily see the smearing related with ray sampling for the broad deepest anomaly. The small heterogeneity is well reconstructed (see horizontal and vertical sections A and C of the figure 15) but the model contains high frequency perturbations. The big one is not well fitted, due to enlightening which is not as dense as in the shallow part. Still results contain a reasonable amount of low frequency information. Nevertheless high frequency information are quite strong and should be removed somehow.

Dealing with the under-determination of the linearized system because of the poor ray coverage could be achieved throuhg smoothing gaussian operators, although these operators are hard to define when working with slowness coefficients. The standard encountered drawback is the over-smoothing by cutting down too small eigenvalues of the linearized system. The level of smoothing needed to remove high frequency perturbation in low resolution zone leads to a loss of high frequency in better defined parts of the model (see Figure 16).

Wavelet decomposition allows different resolutions at different places in the model, as the traveltime tomography predicts. Once the Fréchet derivative is computed, we applied the M operator defined as follows.

From the surface down to the depth of 800 m, wavelet coefficents are kept at all scales. From the depth of 800 m down to the bottom, only the fourth level of wavelet coefficients is kept, corresponding to a resolution of 400 m, coefficients of the finer other three scales are zeroed out. The inversion provides a reconstruction which varies in resolution as we move to deeper depths (Figure 17). We have chosen the mask operator in order to separate zones of the two targets. This illustrates the local resolution with the wavelet decomposition: a high resolution in the shallow part and a broader resolution in depth. One can see than the section A is quite similar to the same section in the figure 15. The vertical section C displays the same profile for the shallow part as the figure 15 while the deepest part is smoother and similar to the same section of the figure 16. Similar features are displayed for the section D: the shallow part still contains high frequency information, as all wavelet coefficients have been kept, whereas the deep part is a smooth version of the same section of the figure 15. The high frequency artefacts have also been removed from the horizontal section B of the largest heterogeneity, thanks to the wavelet decomposition and the mask application.

We may conclude that wavelets enable us to adapt locally the resolution to the theoretical resolution of traveltime tomography. This adaptation is done with a multiplication between wavelet coefficients and the mask values, as the Fréchet derivative matrix is calculated in the wavelet domain. How to choose the mask operator is a problem we face?

5 Resolution analysis for defining the binary M operator

For better focusing of the reconstruction through different scale levels, one has to proceed through a resolution analysis in order to overcome difficulties coming from the linearized inversion. A more exhaustive exploration of the parameter space may be of great interest but it is beyond the scope of this paper. We assume that a model has been designed for which one can perform a resolution analysis. It could be the initial model or any model obtained through previous attempts of inversion. Starting with this model, we may want to have an rough estimation of the resolution matrix for better performance of the inversion [START_REF] Nolet | Explicit, approximate expressions for the resolution and a posteriori covariance of massive tomographic systems[END_REF]. Of course, once done, we may end up with another medium we shall use again for a new resolution analysis and so on. We cannot rule out this loop over adapted models.

For small systems, one can use a direct solver and computing the resolution matrix is quite fast. For important systems, we may rely on numerical iterative techniques as proposed by [START_REF] Latorre | A new seismic tomography of Aigion area (Gulf of Corinth, Greece) from the 1991 data set[END_REF]. As shown in the figure 18, a small local velocity perturbation is introduced with an amplitude of 0.5 km/s, making a spike with a given spatial extension. This amplitude is such that we have a perturbation of travel times above the numerical noise and that rays are kept nearly undeformed. Therefore, we may apply a linear procedure for the reconstruction of this spike. We keep rays and Fréchet derivatives computed in the reference model speeding up dramatically the spike reconstruction as we perform it in different zones of the medium for different spatial extensions of the spike. At a given scale, a resolution matrix estimation as shown in the figure 18 will drive us in the reconstruction of the mask at that scale. We zero out Fréchet derivative when less than 95 % of the wavelet amplitude is recovered at the central position of the spike. This decrease of the amplitude is related to the spreading of the reconstructed spike. Hence, we analyze only the amplitude right at the position where is the spike. Because the system is linear, positive or negative amplitudes give the same resolution matrix extimation.

A realistic crosswell experiment: application to the marmousi model

This procedure has been applied to a synthetic crosswell experiment performed inside the marmousi model [START_REF] Bourgeois | Marmousi, model and data[END_REF][START_REF] Versteeg | The marmousi experience: Velocity model determination on a synthetic complex data set[END_REF]. Between two wells of 3 km depths at a distance of 1 km (Figure 19), ray tracing has been performed for computing synthetic travel times. A smooth model has been constructed

and a resolution analysis performed on this smooth model leading to an estimation of a resolution map with a 95 % threshold level (Figure 19). This resolution analysis has led us to zero out the finest fifth level of Haar coefficients we are using in our wavelet decomposition as well as specific zones of other scale levels (Figure 20). In other words, we may start from a very fine grid and the resolution analysis will tell us if this grid could be used for inversion. Inversion is performed using either slowness coefficients or wavelet coefficients with the mask application (Figure 20). One can see that, where the ray sampling is not fine enough especially at the bottom of the two wells, a coarser description of the reconstructed image has been designed through the resolution analysis. Wavelet procedure seems to perform better and provides nicer features than the slowness reconstruction.

7 Application to a real data set

We have analyzed a real data set with a global offset configuration using a surface-surface geometry. A impressive number of shots as high as 702 shots with a regular spacing of 25 m have been recorded by 1440 receivers along a 18 km 2D profile with a spacing of 12.5 m. We have considering half the sampling over shots. First-arrival travel-times for 301 shots have been picked automatically with an industrial software for the entire receiver line and a tentative macromodel of 18 km long and 4 km deep should be reconstructed using these travel-times. The figure 21 displays one shot gather with initial picked time curve and computed time curve. Our slowness inversion grid has 169 vertical nodes and 725 horizontal nodes with a spatial step of 25 m. We design an initial model for the tomography by fitting a mean traveltime curve with a 1D velocity model. We apply this 1D model below the given topography typical of a foothill structure (Figure 22) using a spatial discretization of 100 m on both directions. The model through FADTT procedure has to fit as much as possible the data and it also has to be smooth especially for deep zones as the ray sampling becomes quite poor.

By a proper gaussian smoothing strategy with depth, a 2D velocity structure has been constructed through slowness coefficients (figure 23) with the same spatial discretization as the initial model. No damping is applied except the internal damping of the LSQR procedure. The number of linearized inversion is five and the rms final error is 26 ms. This model seems to take into account spatial low-frequency contents while the high-frequency ones are missing.

This model is used for a resolution analysis (Figure 18) in order to build up the M mask operator which is displayed in the figure 24. We have selected three levels of the wavelet decomposition starting from 100 m spatial spacing. The coarsest grid is therefore with a spatial discretization of 400 m. Using these three characteristic lengths, we have performed resolution analysis and we have adopted a criterion of a 95 % amplitude reconstruction at the center of the spike for keeping Fréchet derivatives coefficients at that level (Figure 24). One can see that the mask is the identity at the scale of 400 m.

By using linear interpolation wavelet coefficients with two vanishing moments, after thirty iterations, a model with an higher spatial frequency content nearby the free surface has been obtaines (Figure 25). Patterns related with the ray sampling have kept under control and do not overprint the reconstructed image at depth.

The rms final error does not show an significant decrease for assessing the quality of the recovered image through this simple quality criterion. By displaying residual times between sources and receivers (Figure 26), one may notice the difficult fitting at near offsets related both to a difficult picking as different phases cross each other and to the model description which might not be accurate enough for ther local variations often represented by source and receivers statics. During this work, we have not used these quantities known to be quite efficient for describing features under the mesh size [START_REF] Kissling | Initial reference models in local earthquake tomography[END_REF].

One way of assessing the resolution variation could be through a checker-board test applied to the final best resolved model. As shown by the figure 27, different scales are better reconstructed by wavelet decomposition than by slowness decomposition. The checker-board over the smoothest patch with a characteristic length of 2 km provides similar results when considering slowness or wavelet decompositions. Other scales are significantly better recovered using the wavelet decomposition, showing us the interest of this decomposition.

Conclusion

Seismic tomography of first-arrival delayed times could be performed using a wavelet decomposition of the model space. By doing so, one is able to control locally the expected resolution of the inversion. Using the linearized procedure does not automatically inject slowness anomalies at the appropriate scale as the resolution vary during the linearized iterations. Consequently, one must introduce an additional operator in order to prevent spatial aliasing which could lock the inversion procedure. The introduction of an binary operator over scale grids will forbid wavelet coefficients perturbation at a scale when we do not expect resolution. This operator M is a local operator: it will allow us reconstruction of various features as detailed ones nearby the free surface simultaneously with coarse ones at depth for a surface-surface geometry. An illustration on a synthetic crosswell experiment should similar features.

Improvement of data misfits is not significant and may be quite below the noise level which leaves us without clear quality control criterion. We are left with a strategy based on the model space representation and the wavelet decomposition leads to similar results as a local covariance analysis.

Application to a real data set has been performed and shows that one can indeed perform seismic tomography using wavelet formulation over a quite important industrial 2D line.

Of course, the attention is put back into picking strategy as we have found that industrial softwares might misdetect picked times and one has to find alternative strategies for better picking as we have now accurate tools for better reconstruction of detailed features through travel-time tomography.

Extension to 3D configuration should be performed and combining this wavelet decomposition and the finite frequency tomography [START_REF] Montelli | Finite-frequency tomography reveals a variety of plumes in the mantle[END_REF] will certainly improve in the future local resolution of seismic imaging before tackling amplitude analysis of seismic traces [START_REF] Pratt | Gauss-newton and full newton methods in frequency-space seismic waveform inversion[END_REF][START_REF] Operto | Quantitative imaging of complex structures from multi-fold wide aperture seismic data[END_REF]. V

1 V 2 V 0 V 3 V 5 V 4

Figure 1 :Figure 2 :

 12 Figure 1: Projection of a 1D signal over 3 levels of wavelets

Figure 3 :

 3 Figure 3: Decomposition of a V 0 slowness model for a crosswell experiment onto the V i spaces using 5 levels of wavelets

Figure 5 :Figure 6 :Figure 7 :Figure 8 :

 5678 Figure 5: Synthetic positive anomaly to be reconstructed: the perturbation model and the ray coverage

Figure 9 :

 9 Figure 9: Synthetic example: a circular negative anomaly on top of an homogeneous model. Please, note the different ray coverage as for a positive anomaly essentially in the center of the box. The value should not be recovered in the center and only the parametrization makes this part of the model sensitive to travel times.

Figure 10 :

 10 Figure 10: Tomographic inversion without any regularization or damping. On the top left panel, the reconstruction anomaly and the ray sampling. On the top right panel, the misfit behaviour leading to a minimal rms error of 1.50 ms. Four iterations are enough for convergence. Please note high frequency noise and the glitch along geometrical axes related to poor ray coverage as well as edge effects for similar reasons.

Figure 11 :Figure 12 :Figure 13 :

 111213 Figure 11: Tomographic inversion using wavelet with two vanishing moments and two levels. The convergence is obtained with an rms error of 1.24 ms after ten iterations. See figure 10 for other comments

Figure 14 :

 14 Figure 14: The synthetic model with a gradient in velocity and two circular anomalies on top panel (one is positive and the other negative). The model panel shows only these two anomalies with a central amplitude of 0.5 km/s while the bottom panel displays the ray coverage in the true synthetic model.

Figure 15 :Figure 16 :

 1516 Figure 15: Reconstruction in the slowness domain whithout any regularization. Please note the ray sampling imprint into the reconstructed image as we want high resolution for the reconstruction of the small superficial anomaly

Figure 17 :

 17 Figure17: Reconstruction using wavelets with four vanishing moments and four levels: a mask ensures that, under the depth of 800 m, the fourth level of wavelet is the only one kept corresponding to a resolution of 400 m. The large heterogeneity is kept smooth and the shallow heterogeneity is well reconstructed.

Figure 18 :

 18 Figure 18: Resolution analysis: top panel displays the spike anomaly we want to reconstruct, the medium panel the spreading of the amplitude reducing it at the center and the bottom panel the mask operator as the percentage of the recovered amplitude.

Figure 19 :

 19 Figure 19: Resolution analysis: both two panels on the left display the true model between two hypothetical wells sampling this zone of the Marmousi model with quite dense ray coverage, the third panel shows the initial smooth model on which the resolution analysis is performed while the fourth panel shows resolution map as percentage of the recovered spike amplitude.

Figure 20 :

 20 Figure 20: Tomographic inversion: the left panel displays the mask operator over the five scale levels from the finest one equivalent to the sampling of the synthetic Marmousi model to the coarsest one, the middle panel the inverted model using slowness coefficients without any damping and the right panel the inverted model using Haar wavelet coefficients over five scales.Please, note that cancelling the entire finest scale level does not harm the inversion procedure.

Figure 21 :

 21 Figure 21: First-arrival picked time using an industrial software for the gather of shot 311; synthetic travel-times are displayed as well showing capacities of interpretation through the recovered velocity model

Figure 22 :

 22 Figure 22: Initial 1D model for the travel-time inversion: this model starts exactly at the free surface following a complex topography.

Figure 23 :

 23 Figure 23: Reconstruction through slowness coefficients inversion: smooth gaussian regularization has been applied, preventing details to be exhibited as well as ray smearing features.

Figure 24 :

 24 Figure 24: Mask construction over the already constructed smooth medium. Please note that the finest scales are only active nearby the free surface

Figure 26 :Figure 27 :

 2627 Figure 26: Final time residuals after thirty iterations. Please note the still too fast velocity structure nearby the free surface.Very local description of a low-velocity zone should be given and could be performed through statics evaluation.

km, the left panel is for slowness coefficients while the right panel is for wavelet coefficients. Please note the better behaviour through wavelet coefficients inversion.

Acknowledgments: This work is done under the industry French ministry contract FSH N • CEP&M RE.1006/05. We thank TOTAL for providing us the complete data set as well as the industrial picked time.

Many interesting discussions with Frank Adler of TOTAL have allowed the improvement of this work.