Philippe Maurine
email: pmaurine@lirmm.fr

Nadine Azemard
email: azemard@lirmm.fr

Daniel Auvergne
email: auvergne@lirmm.fr

N Azémard

Deep Submicron Switching Current Modeling for CMOS Logic Output Transition Time Determination

Non zero signal rise and fall times contribute significantly to CMOS gate performances such as propagation delay or short circuit power dissipation. We present a closed form expression to model output rise and fall times in deep submicron CMOS structures. The model is first developed for inverters considering fast and slow input ramp conditions. It is then extended to gates through a reduction procedure considering the maximum current available in the serial transistor array. Validation of this modeling is obtained by comparing calculated gate output transition time to simulated ones (HSPICE level and foundry card model on 0.18µm process).

INTRODUCTION

The use of safe gate level characterization of performances over the full design space is the only way to maintain timing relationships between functional blocks when designs approach complexity of millions of transistors. To control or drive design alternatives, technology migration, as well as process variation it appears necessary to get available design oriented models to evaluate the performances of specific structures. The traditional representation of delay associates a constant "inertial" delay characteristic of the cell to an output load dependent delay characterizing the cell size and structure.

However input-to-output coupling effects associated to speed saturation of the carriers induce non linearity for the propagation delays which are important enough to be considered for accurate cell delay-performance characterization. Great sensitivity of the delay to the edge of the input controlling signal has been observed in submicron processes. These edges are generally defined as the controlling gate output -voltage transition time measured between appropriate voltage levels. These signal rise and fall times contribute significantly to the delay and are responsible of the nonlinear variation of real delay values. As a result, gate delay characterization implies consideration of propagation and output transition times.

The modeling of the gate output transition time has been the object of numerous works. Due to the difficulty in solving the complete differential equation representing the discharge (charge) of the gate output node, various attempts have been done to characterize this output transition time, including step [START_REF] Burns | Switching response of complementary symmetry MOS transistor logic circuits[END_REF], ramp [START_REF] Hedenstierna | CMOS circuit speed and buffer optimization[END_REF] and exponential models [START_REF] Kayssi Ayman | Analytical transient response of CMOS inverters[END_REF]. In [START_REF] Sakurai | Alpha-power model, and its application to CMOS inverter delay and other formulas[END_REF], a submicron delay and output slope modeling is given, still limited to fast input transitions. Recently, as an extension of the work proposed in [START_REF] Sakurai | Alpha-power model, and its application to CMOS inverter delay and other formulas[END_REF], S. Dutta [START_REF] Dutta | A Comprehensive Delay Model for CMOS Inverters[END_REF], considered very slow input ramp effects. Both the delay and the output ramp duration are obtained by curve fitting between two extreme points corresponding to infinitely fast and infinitely slow inputs. As an improvement of his initial Deep Submicron Switching Current Modelling for CMOS Logic 5.3.2 work Sakurai [START_REF] Sakurai | A simple MOSFET model for circuit analysis[END_REF] considered extremely fast and slow ramp conditions and solved intermediate cases from smooth interpolation between the two extremes. In [START_REF] Bisdounis | Analytical transient response of propagation delay evaluation of the CMOS inverter for short channel devices[END_REF] Bisdounis proposed a fast and slow input slope definition from the operating mode of the switching transistor however, no clear design oriented definition of both fast and slow input transition range, based on the size and the load of the switching and controlling devices appears available. Hirata in [START_REF] Hirata | Proposal of a Timing Model for CMOS Logic Gates Driving a CRC π load[END_REF] proposed a piece wise linear representation of the current available in the switching structure. This approach necessitates a great number of calibrations with Spice simulations of the different technological parameters used in the representation. In fact the output ramp duration of a CMOS structure depends on its current possibility (I MAX ) and of the amount of charge to be transferred (C.V DD ). As proposed in [START_REF] Auvergne | Signal transition time effect on CMOS delay evaluation[END_REF] it can be obtained from:

MAX DD L OUT I V C . = τ ( 1
)
where V DD represents the node voltage variation and C L it's output loading capacitance.

As shown the key parameter in modeling the output transition time is the current available in the switching structure of which determination depends on the structure, its size and the duration time of the input controlling edge. In order to complete an analytical model of delays developed for submicron CMOS structures [START_REF] Auvergne | Signal transition time effect on CMOS delay evaluation[END_REF], we present in this paper a design oriented macro modeling of the CMOS structure output transition time. In section 2 we present the method we used to obtain the value of the maximum current available in CMOS inverter and gates considering both fast and slow input ramp conditions. The modeling and the validation of output transition time is given in section 3. Section 4 draws a conclusion on this model.

INVERTER MAXIMUM CURRENT

Depending on the strength of the controlling structure two design conditions have to be considered, fast and slow input ramp conditions. Let us consider an inverter with a load C L controlled by a rising linear input ramp of duration τ IN . As shown in Fig. 1, the current sunk from the load by the N transistor depends on the value of τ IN :

-in region 1 the set up of the current of the N transistor follows the input ramp variation and exhibits a constant maximum value during all the discharge process, this defines the fast input range, -in region 2 the maximum current is obtained before the input ramp reaches its maximum value, resulting in a smaller value of the charge evacuated by unit time. This defines the slow input range where the maximum value of the discharging current decreases when the input transition time increases.

Maximum current value for fast input range

During all the input ramping process the N transistor is saturated, its current maximum value is defined for V IN = V DD , resulting in: ) .( .

TN DD N N fast MAX V V W K I - = (2)
where K N is the transistor conduction factor defined in [START_REF] Sakurai | Alpha-power model, and its application to CMOS inverter delay and other formulas[END_REF] for α=1, V TN and W N the N transistor threshold voltage and width respectively.

Maximum current value for slow input range

As in the preceding case the transistor is still in saturation when its current reaches the maximum value but its gate driving voltage is smaller and its value must be defined. For that we consider that in the time interval t VTN -t MAX , (Fig. 2), the current exhibits a linear variation. This gives:

) . .( . ) (

TN IN DD N N N V t V W K t I - = τ (3)
and:

IN DD N N MAX V W K t I t I τ . . = ∆ = ∆ ∆ ( 4
)
where the input ramp duration time τ IN is the output ramp duration of the controlling structure, as defined in eq.1.

Under the approximation that the current variation is symmetric with respect to its maximum value we can evaluate the total charge removed at the output node as:

2 . 2 . t I V C MAX DD ∆ = (5)
Combining eq. 4 and 5 we obtain the value of the maximum current resulting from a slow rising input controlling edge as:

IN DD N N slow MAX C V W K I τ . . . 2 = (6)
where C.V DD represents the total charge to be removed from the output node, where: C=C L +C SC +C PAR in which C L and C PAR represent the inverter active load (output loading gates) and the output parasitic capacitance respectively, C SC is the short circuit equivalent capacitance which represents the charge by volt unit between the supply rail during the discharge process as defined in [START_REF] Daga | A comprehensive delay macromodeling for submicron CMOS logics[END_REF][START_REF] Turgis | A novel macromodel for power estimation for CMOS structures[END_REF]. We compare in Fig. 3 the maximum current values deduced from eq.2 and 6 to the values simulated with Hspice for an inverter defined by W N =1µm, W P =2.2µm, L=0.18µm for different loading conditions (5,10,15 and 20 times its input capacitance C IN =4.5fF). As shown we obtain a very good agreement between simulated and calculated values (less than 10% discrepancy) over the considered full design range.

2-3 Maximum current value for a simple gate

To evaluate the maximum current available in a gate it is necessary to consider the current limitation effect produced by the serial array of transistors, together with the multiplication effect produced by the dual parallel array. The current possibility of this parallel array is input vector dependent, but bounds can be easily defined considering, for an n input gate, one or n times the maximum current of an inverter with identically sized transistors. The reduction of the serial array to an equivalent transistor has been the object of numerous works [START_REF] Hirata | Proposal of a Timing Model for CMOS Logic Gates Driving a CRC π load[END_REF][START_REF] Chatzigeorgiou | Collapsing the Transistor Chain to an Effective Single Transistor[END_REF][START_REF] Nabavi-Lishi | Inverter Models of CMOS Gates for Supply Current and Delay Evaluation[END_REF][START_REF] Wang | A New Short Circuit Power Model for Complex CMOS Gates[END_REF]. To reduce a gate to an equivalent inverter we present here a new reduction method by considering the serial array of n transistors as an input voltage controlled current generator, as illustrated in Fig. 4.

If we consider a control on the top input (Bot and Mid inputs connected to V DD ) of the Nand3, we can see easily that the voltage dropt through the Mid and Bot transistors, working in linear mode, reduces the voltage swing of the controlling gate. This results in a transistor size dependent reduction of the available current in the network with can be modeled as a reduction factor equal to the ratio of the currents available in the array and in the inverter with identically sized transistors. This gives:

N N N fast .R .W K d + =1 Re ( 7
)
where R N represents the sum of the resistance of the bottom transistors. For a control on the bottom input, for fast input ramps, the intermediate nodes are discharged faster than the output one. In this case the current is still limited by the top transistor and the reduction factor is given by eq.7.

For slow input ramp condition the bottom and top transistors operate in saturated mode and the current is limited by the bottom transistor working with a reduced drive and drain source voltage. In this condition it appears necessary to calibrate, from simulations on the process, the conduction factor of the bottom transistor in the serial array [START_REF] Nose | Analysis and future trend of short circuit power[END_REF]. For the process under study (0.18µm) values of Red SLOW = 1.2, 1.48 and 1.78 have been obtained for NAND 2, 3 and 4 respectively, which are quite different from the values obtained for fast edge conditions (1.55, 2.1 and 2.6 for NAND 2,3,4 respectively) or from a direct reduction based on the number of serial transistors [START_REF] Sakurai | Alpha-power model, and its application to CMOS inverter delay and other formulas[END_REF].

Controlling the middle input we obtain a superposition of the contributions of the preceding effects. The reduction factor can easily be deduced from the preceding cases considering the middle transistor in top or bottom position for the bottom or top transistor of the array, respectively, resulting in a reduction factor:

SLOW FAST d d d Re . Re Re = (8)

OUTPUT TRANSITION TIME

The output transition time can be obtained easily from eq.1 by replacing I MAX by the expressions previously developed.

Inverters

Considering fast and slow input ramp conditions results in:

          - = fast OUT IN DD TN DD fast OUT OUT V V V MAX τ τ τ τ .
; [START_REF] Auvergne | Signal transition time effect on CMOS delay evaluation[END_REF] with:

HLS N L ST fast OUT T C C . 2 . = =τ τ (10) ( ) N TN DD OX GEO DD ST K V V C L V . . . - = τ ( 11
)
In these equations T represents the step response of the inverter, and τ ST the shorter switching time of the process, as defined in [START_REF] Daga | A comprehensive delay macromodeling for submicron CMOS logics[END_REF]. Validation of these expressions has been realized on different configurations of inverters in various loading and controlling conditions by comparing simulated (Hspice BSIM3 level 49) and calculated (eq.9) output duration time values. The results obtained are illustrated in Fig. 67. The output transition time evolution is given versus the ratio τ IN /T HLS used as a metric for input transition times. The expression for an output rising edge can be obtained by exchanging N and P suffixes. As shown we obtain a very good agreement between simulated and calculated values (less than 10% discrepancy) over the considered full design range.

Gates

Considering the current reduction factors defined in eq.7-8, the generalization to gates is straightforward, we obtain:

for a Top input control:

( )           - = fast OUT IN DD TN DD FAST fast OUT FAST OUT V V V d d MAX τ τ τ τ . . Re ;
. Re [START_REF] Chatzigeorgiou | Collapsing the Transistor Chain to an Effective Single Transistor[END_REF] for a Bot input control:

( )           - = fast OUT IN DD TN DD SLOW fast OUT FAST OUT V V V d d MAX τ τ τ τ . . . Re ;
. Re [START_REF] Nabavi-Lishi | Inverter Models of CMOS Gates for Supply Current and Delay Evaluation[END_REF] for a Mid input control:

( )           - = fast OUT IN DD TN DD fast OUT FAST OUT V V V d d MAX τ τ τ τ . . . Re ;
. Re [START_REF] Wang | A New Short Circuit Power Model for Complex CMOS Gates[END_REF] Validation has been done following the same procedure than for inverters.

Discussion on slow and fast ranges of input transition times

Valuable criteria in evaluating the quality of designs or in defining metric for design performance optimization is to clearly identify the limit condition between fast and slow transition times. As previously illustrated this limit (fig. 1) corresponds to the threshold between the availability of constant discharging (charging) current and varying one. It depends of the relative values of the input and output transition times.

For example let us identify this limit for inverters in equalizing the two terms of eq.9, (eq.12-14 for gates). This gives the limit at which input ramps must be considered as slow as:

FAST TN DD DD IN OUT V V V τ τ .         - ≥ (15)
Remembering that on an array τ IN represents the input transition time of the controlling inverter (i-1) and τ OUT the output transition time of the switching device (i) we obtain from eq.15:

) ( ). ( ) 1 ( i Fo R k V V V i Fo DD TN DD µ - ≥ - (16)
) ( ). ( ) 1 ( i Fo k R V V V i Fo DD TP DD µ - ≥ - (17)
for output falling and rising edges, respectively, where Rµ is the ratio of the conduction factors of N and P transistors, k and Fo have been previously defined.

Extension to gates can be easily obtained from eq.12-14, including the reduction factors.

In table 2, we compare the limit value of the τ IN separating fast and slow input range as defined in eq.15, to the values deduced from the simulation on an inverter for different configuration ratio values and loading conditions. As shown the limit previously defined is in very good agreement with the values obtained from the simulations. In the figure 8 we illustrate the relative character of the definition of the limit between fast and slow input ranges. The curves represent the output voltage and discharging current controlled by an input ramp of 50 and 500ps of duration with different output loads. As shown, (fig. 8.a) while the input duration ramp is quite short (50ps), due to the weakness of the load (Fo=0.5) the input control must be considered as slow. On the other hand (fig8.b) a heavy loaded inverter (Fo=80) controlled by a quite long duration ramp (τ IN =500ps) is controlled under fast input ramp conditions. This results is very important, this justifies why, as well for defining design validation range than look up tables, it is necessary to define the input control range relatively to the output transition time of the considered cell.

CONCLUSION

We derived design oriented simple and closed form formula for the output transition time of CMOS gates. We showed that the proposed expression reproduces the sensitivity to the design and process parameters. Based on a metric defined on inverter for fast input ramp conditions the formula includes deep submicron effects by considering the variation of the maximum current available with the input edge. Extension has been done to gates by reduction to an equivalent inverter, considering the different input control conditions. Clear evidence of different reduction factor values for fast and slow input edges has been demonstrated. Validations through Hspice simulations for a 0.18µm process confirmed the validity of the proposed expressions which can easily be used to replace look up tables in timing estimator.

Clear definition of the slow and fast input control range is clearly defined and demonstrated. Application to edge control for low power buffer design is under development.

Fig. 1 .

 1 Fig. 1. Illustration of the fast ¬ and slowinput controlling ranges of an inverter.

Fig. 2 .

 2 Fig. 2. Discharging current evolution. Fig. 3.Comparison between calculated and simulated maximum discharging current value; ¬ andlabel the fast and slow input ranges.

Fig. 4 .

 4 Fig. 4. Reduction of the serial array of transistors to a multiple input voltage controlled current generator.

Fig. 5 .

 5 Fig. 5. Static I/V characteristics of the current generator with respect to the controlling input (Top, Mid, Bot).

Fig. 6 .Fig. 7 .

 67 Fig. 6. Inverter output transition time (W N =1µm, k=2, L=0.18µm) loaded by 5 to 20 C IN .

Fig. 8 .

 8 Fig. 8. Illustration of the relative definition of slow ((a) τ IN =50ps) and fast ((b) τ IN =500ps) input range, obtained with short and long duration ramps respectively.

HLS SIM CAL ∆% SIM CAL ∆% SIM CAL ∆% SIM CAL ∆%

  Table 1 and 3 are relative to Top and Bottom controlled Nand2,3 (W N =W P =1µm) loaded by 10.C IN and implemented in a 0.18µm process. As shown we obtain a very good agreement between simulated and calculated values of the output transition time.

	τ IN /T 2	Nand2_Top 144 154 7% 240 223 5% 134 121 10% 197 204 6% Nand3_Top Nand2_Bot Nand3_Bot
	6	178 172 3% 256 252 8% 166 155 7% 209 204 0%
	10	227 222 2% 315 325 3% 201 199 1% 233 217 1%
	16	291 281 3% 392 385 2% 275 268 3% 296 291 5%
	20	329 315 4% 439 436 2% 292 282 4% 312 306 4%

Table 1 .

 1 Comparison between simulated and calculated values of output transition time for NAND2, 3 with top and bot input control.

	Deep Submicron Switching Current Modelling for CMOS Logic	5.3.8

Table 2 .

 2 Comparison between simulated and calculated (eq.15) values of the limit value between fast and slow range defined for τ IN .

	τIN limit	Fo=5		Fo=20
	(ps)	Sim	Cal	Sim	Cal
	k=1	63.1	60	243	240
	k=2	94.6	90	365	360
	k=3	126	120	487	480


Mission Lu Pan

Cathy Quantin

Chloé Michaut

THE COMPOSITION AND STRATIGRAPHY OF THE NORTHERN LOWLANDS AND IMPLICATIONS FOR INSIGHT

Introduction:

The Martian dichotomy presents a striking contrast between the northern and southern hemispheres in elevation [1], crustal thickness [2], surface roughness [3] and composition [4,5]. Data acquired in the recent decades give clues to diverse aqueous environments [4,5] and igneous lithologies [6,[START_REF] Cousin | Icarus[END_REF] in the southern highlands, while much is left to explore within the northern lowlands. The subsurface of the lowlands records the major stages in Mars evolution, including the formation of the Martian dichotomy, the evolution of the Martian crust, the existence of a global ocean, etc. Here we present the updated understanding of the geologic processes in the lowlands and discuss its implications for the InSight (Interior Exploration using Seismic Investigations, Geodesy and Heat Transport) mission to Mars.

A hemispheric dichotomy of the Martian crust?: Crustal thicknesses in the highlands and lowlands interpreted from topography and gravity models [1,2] rely on the assumption of uniform bulk density. However, if the northern lowlands formed from a giant impact, the lowlands crust could be partially formed from subsequent mantle melts, with a distinct composition from the highlands. Based on constraints from the Martian elastic lithosphere thickness, thermal evolution models suggest the highlands could have higher radiogenic elements concentration and smaller bulk density than the lowlands [START_REF] Thiriet | [END_REF], which would result in different crustal thickness and heat flow estimates [START_REF] Thiriet | [END_REF][START_REF] Wieczorek | Mars crust conf[END_REF].

The largest craters may excavate materials beneath the sedimentary layers in the northern lowlands and expose part of the underlying crust. Through investigations of impact craters using data acquired by CRISM (Compact Reconnaissance Imaging Spectrometer for Mars), diverse hydrated minerals (including Fe/Mg phyllosilicates, chlorites, silica) have been identified, indicating a hydrated basement overlain by mafic mineral bearing units [START_REF] Carter | [END_REF]11]. The mineralogical classes imply a similar, altered basement as the southern highlands, and there is yet no strong geological evidence for different crustal materials.

Mineralogical record of sedimentary and volcanic infilling: Earlier study using data acquired from Thermal Emission Spectrometer (TES) showed the lowlands surface differ from the basaltic type terrain in the highlands [12] and the spectral features are later interpreted to be due to a veneer of alteration products [13,14] on the surface of the lowlands. Through the mineralogy exposed in impact craters, widespread mafic minerals have been identified [11,15], confirming the hypothesis that the lowlands are filled with wrinkle-ridge forming mafic lava flows [START_REF] Head | JGR Planets[END_REF]. The identification of layered deposits in Chryse and Acidalia Planitia [START_REF]Geology[END_REF] supports the notion that sedimentation occurred after the volcanic emplacement, possibly due to late outflow channel activities. Interestingly, no direct evidence of widespread carbonates or chlorides has been identified in the lowlands [11].

Implications for InSight: The complex resurfacing stratigraphy at the dichotomy boundary will pose complexities on the signals of the seismic waves collected by the SEIS (Seismic Experiment for Interior Structures) instrument on board InSight lander. Characterization of the geological context is thus important to the data interpretation. From recent mapping work [START_REF] Pan | [END_REF], we have identified mafic mineral detections surrounding the landing site ellipse, confirming the composition of the volcanic plains unit. Layered sediments found in several crater central peaks suggest a possible extension of the sedimentary unit (Medusae Fossae Formation) to the north and east of the landing site. On the other hand, with the advent of InSight mission to Mars, we may be able to detect the seismic discontinuities in the subsurface of the landing site. The future measurements of heat flow and seismic properties will place major constraints on the thermal and density structure of the crust and shed new lights on the nature of the dichotomy boundary and the complex stratigraphy within the northern lowlands of Mars.

Acknowledgement: This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 751164.


Smooth velocity models in reflection tomography for

imaging complex geological structures.


Matthieu Armando

Jean-Sebastien Franco

Edmond Boyer

Adaptive Mesh Texture for Multi-View Appearance Modeling

In this paper we report on the representation of appearance information in the context of 3D multi-view shape modeling. Most applications in image based 3D modeling resort to texture maps, a 2D mapping of shape color information into image files. Despite their unquestionable merits, in particular the ability to apply standard image tools, including compression, image textures still suffer from limitations that result from the 2D mapping of information that originally belongs to a 3D structure. This is especially true with 2D texture atlases, a generic 2D mapping for 3D mesh models that introduces discontinuities in the texture space and plagues many 3D appearance algorithms. Moreover, the per-triangle texel density of 2D image textures cannot be individually adjusted to the corresponding pixel observation density without a global change in the atlas mapping function. To address these issues, we propose a new appearance representation for image-based 3D shape modeling, which stores appearance information directly on 3D meshes, rather than a texture atlas. We show this representation to allow for input-adaptive sampling and compression support. Our experiments demonstrate that it outperforms traditional image textures, in multi-view reconstruction contexts, with better visual quality and memory footprint, which makes it a suitable tool when dealing with large amounts of data as with dynamic scene 3D models.

Introduction

Image based 3D shape modeling is the process of building digital models of shapes using real images. It finds applications in many domains, in particular with the new virtual and augmented reality devices and the associated need for 3D contents. In order to represent the reconstructed 3D shapes, the dominant paradigm is to model them as geometric surfaces over which appearance functions are defined. With such a model, both geometric and appearance features fundamentally contribute to convey realism and fidelity to † MSR-INRIA Joint Center Figure 1: Image-based 3D modeling: given a set of input photographs (left), a geometric mesh is computed (top), along with an appearance function stored within the surface mesh structure (bottom), the main focus of this paper. the observed shapes. In this work we focus on the appearance representation, which has received surprisingly little attention from the 3D vision community. In particular, we show that both visual quality and encoding efficiency can be significantly improved by exploring representations beyond traditional 2D texture maps.

Assume we are given a geometric shape model estimated from n images, using for instance a multi-view stereo approach [START_REF] Leroy | Shape reconstruction using volume sweeping and learned photoconsistency[END_REF]. An optimal solution for the appearance, and with respect to the observed information, is to keep all the original images in the representation. This appears to be inefficient, or even intractable, in many situations where possibly numerous high resolution cameras, and potentially temporal sequences of shapes, are considered. Hence, most 3D visual modeling frameworks usually resort to an appearance representation in the form of a 2D regular grid of texels whose values, typically RGB colors, are estimated from the observed images, given the geometry. The main benefits of these texture maps are their controlled and limited sizes in addition to inheriting most regular 2D image tools, including compression and filtering. On the other hand, such representations present severe limitations. First, a 3D to 2D mapping is required to associate geometric points on the 3D shape to appearance information in the 2D texture. To this purpose, the 3D shape model is usually cut into charts homeomorphic to discs that are parameterized with texel coordinates. The resulting 2D texture atlases are discontinuous by construction, creating unnecessary seams on the model appearance and making local computations on the appearance function difficult; In addition, locally adjusting the per-triangle texel sampling to the density of pixels oberving each triangle is unpractical, as it is globally bound by the initially chosen texture resolution and would require global optimization of the atlas mapping function.

Instead of an appearance stored in an intermediate 2D grid, we advocate therefore for an appearance stored on the 3D shape surface with the following advantages:

• There is no need for a 3D to 2D mapping that induces discontinuities in the representation.

• Regions on the shape with varying observed pixel densities, as due to camera properties including resolutions and viewpoints, as well as to self-occlusions such as armpits on the human body, can be described with adapted appearance samplings.

To this end we build on the mesh color representation introduced in [24] for synthetic graphical models and propose an adaptive appearance representation for 3D image based reconstruction. We show that this representation outperforms traditional texture atlases with more efficient appearance representations with respect to both precision and size. In addition, we also introduce a dedicated compression method that builds on the JPEG pipeline and allows for effective appearance storage with our representation. The efficiency improvements of the proposed contributions allow for promising image based 3D modeling applications.

Related Work

Appearance Representation.

Besides the widespread texture map representation mentioned before, attempts have been made for representations beyond 2D textures in order to address their limitations. The thorough survey of Tarini et al. [START_REF] Tarini | Rethinking Texture Mapping[END_REF] gives an exhaustive overview on the subject. For instance Ptex [3] used by Walt Disney Animations Studios, and Mesh Colors [24] are particularly noteworthy as they are used in real production pipelines in the movie industry ([10]). In both cases, they use the intrinsic parameterisation of the model, rather than an arbitrary mapping. While Ptex was developed for subdivision surfaces, and only works with quad faces, Mesh colors can be used for triangular meshes which is the 3D model representation we consider in this work. However, works on this model focus yet on the graphics side, either trying to optimize it for rendering [23] or creating 3D painting tools for artists to work with [10]. For image-based modeling in vision applications, texture maps are still the main appearance representation and, to the best of our knowledge, no attempts have been made to replace them.

Multi-view Appearance Estimation.

Appearance representations are used to store the appearance information for 3D models. In the multi-view reconstruction context, such information is estimated by considering the color cues from several input views. In the case of texture maps, texels are usually backprojected into input images and the corresponding pixel values are combined with: blending approaches that average pixel colors to determine texel values, e.g. [16,[START_REF] Baumberg | Blending Images for Texturing 3D Models[END_REF]; or more recently using superresolution strategies, as in [9, [START_REF] Tsiminaki | High Resolution 3D Shape Texture from Multiple Videos[END_REF]8], that allow for prior constraints, e.g. TV regularization, to be applied on the estimated texture. In any case, the disjoint charts of a texture atlas introduce discontinuities in the texture space that make filtering, regularization, or any non-local computation over the texels difficult to perform, since neighbors must be determined across chart boundaries. In that respect, a representation attached to the shape parametrisation, e.g. mesh triangles, is much more efficient with connectivity relationships that are intrinsic. In addition, it enables for adaptive resolutions where a texture has a uniform pixel resolution.

Appearance Compression.

Reducing the memory footprint is an important issue when streaming or storing 3D contents. This is even more critical with dynamic scenes for which shape and appearance information evolve over time. In the case of 2D textures, appearance information can be optimally compressed using image techniques with nevertheless some differences. Whereas most image compression techniques are optimized for storage and transmission, and thus, focus on visual quality vs. bitrate, texture compression puts more emphasis on random access and decoding speed [2]. State of the art methods include ASTC (adaptive scalable texture compression) [14] and ETC (Ericsson Texture Compression), published as iPACKMAN [START_REF] Ström | i PACKMAN: high-quality, low-complexity texture compression for mobile phones[END_REF]. Putting aside the concern about random access and decoding speed, and focusing on bit rate instead, state of the art methods are then mostly standard image compression techniques such as JPEG. This is extended to dynamic scenes in [4] where the standard H.264 compression format is used to compress texture frames over time, assuming for that purpose the texture atlas to be fixed over frames. This enables temporal redundancy to be exploited, which is an interesting feature however beyond the scope of this paper. With mesh colors representations, for which appearance information does not live in 2D images but on the model representation, i.e. mesh triangles, standard image compression techniques do not directly apply. We propose a strategy that builds on JPEG to compress mesh color representations and to benefit therefore from its visual quality properties without sacrificing memory space with respect to 2D textures.

Appearance Models

Without loss of generality we assume that shapes are modeled with 3D triangle meshes. As mentioned earlier their appearances, typically color information, can be represented with 2D image textures, the traditional model, or directly on shapes with mesh textures. In this section we present both models, their parameterizations and sampling properties before considering more specifically in the next section, the 3D image based modeling context and the mesh texture tools we introduce for that purpose.

Notation

We consider the shape surface M, as obtained with for instance a multi-view stereo reconstruction approach. M is represented as an ordered set of vertices {v 1 . . . v N } in R 3 , and an ordered set of triangles {T 1 . . . T M }, where T i ∈ [1, N ] 3 is a triplet of vertex indices. We want to represent the appearance of M as a dense function C : M → C that associates a color from a space C to every point on the mesh. Typically, C = [0 . . . 255] 3 for a discrete RGB-space color representation. To this aim we parameterize point locations on M with barycentric coordinates and through the following function F :

F : {1, ..., M }×]0, 1[ 2 → M (j, a, b) → a v x + b v y + (1 -a -b) v z , (1)
where (x, y, z) = T j are the 3 vertex indices of triangle T j . The closure of the codomain of F is M. We call F -1 the inner parameterization of M. In other words, every point on M can be uniquely represented by its triangle index j and its barycentric coordinates (a, b), with the exception of vertices and edges. The function C : {1, ..., M }×]0, 1[ 2 → C associates then a color to any point represented by (j, a, b).

Image Texture

We recall in this section the widely used image texture solution and discuss its properties.

Parameterization.

Image textures make use of an external mapping between an image and the 3D mesh M under consideration. M is divided into k patches P 1 , ..., P k ⊂ M, such that ∪ k j=1 P j = M. For each patch, a separate function φ i : 3). The chart building process is thus typically a non trivial and global task, introducing discontinuities such that the whole mapping M → T is piecewise continuous. The mapping φ = ∪ k j=1 φ j is typically defined per face, using 3 pairs of texture coordinates in T representing the images of its 3 vertices by φ. φ is then linearly interpolated inside the face. Vertices and edges are duplicated along chart discontinuities. In other words, a discrete function φ 0 : {v 1 . . . v N } → T is defined on the vertices, and is then extended to the whole surface using the inner parameterization F -1 (1). An additional drawback of the scheme is thus that it must ensure that color information C, or any other surface function defined over T , should be consistent across charts to avoid visual artifacts and apparent color seams.

P i → C i ⊂ R 2 is computed independently that maps 3D points in P i to 2D points in a chart C i . Charts are chosen so that ∪ k j=1 C j = T ⊂ [0, 1] 2 and C i ∩ C j = ∅, ∀i = j (see Fig.

Sampling.

T is discretized with a rectangular image grid, inheriting all the associated image tools. The sampling depends on vertical and horizontal resolutions and the intially chosen mapping of triangles to texture space. Some texels (i.e. pixels in the texture) inevitably fall partly within and outside the chart union T . This means C must be extrapolated over the chart borders to avoid artifacts. This is even more critical if consistent filtering operations on C are to be performed.

Storage and Compression.

Image textures can be stored as regular image files and buffers, with numerous possible compression schemes, such as JPEG. The mapping function φ 0 is usually stored into the geometry file, in the form of texture coordinates. Assume E, N, M are respectivley the number of edges, vertices and faces of M. As shown in [START_REF] Peng | Technologies for 3D mesh compression: A survey[END_REF], for large triangular meshes, the approximation M ≈ 2N is considered valid.

In general, the geometry information is represented by N 3D points, or 3N floating point numbers, and the connectivity is stored as M triplets of vertex indices, i.e. ≈ 6N integers. For appearance, 3M texture coordinate pairs are used, i.e. 12N floating points in the worst case scenario, or four times the geometry information, although this can be optimized with smart indexing. In the best case, neglecting cross-boundary chart redundancies, this takes 2N floating points (one pair of texture coordinates per vertex).

Mesh Texture

In contrast to image textures, the appearance information is stored on the mesh structure with mesh textures. Consequently, a 2D mapping function is not required and the appearance information is directly sampled over the mesh triangles, as illustrated in Figure 4. These ideas were introduced in [24] and we experiment and extend them with appearances from real images.

Parameterization.

Mesh textures do not require any form of external parameterization. Color information is defined directly on the mesh structure and not in any other domain. As such, it only makes use of the inner parameterization F -1 of M (1).

Sampling.

Sampling is a strong feature of the model governed by a single resolution parameter. Given a triangular face T i and its resolution R i , a mesh texture stores (Ri+1)(Ri+2) 2 samples on the locations given by:

F (T i , m R i , n R i ) | 0 ≤ m ≤ R i , 0 ≤ n ≤ R i -m . (2)
In other words, the positions of appearance samples within a triangle are parametrized by the barycentric coordinates

( m R i , n R i , 1 - m + n R i ) , (3)
associated to the triangle vertices. In practice, vertices and edge samples are shared between adjacent faces, and they must be treated separately. The resolution of an edge is defined as the lowest resolution of its two adjacent faces. As in [24], face resolutions can only take values that are powers of two. This makes interpolation easier along edges with a lower resolution than the face. This choice of representation has several important advantages over image texture sampling. First, the sampling is by construction hexagonal, and hence more tightly packed and less directionally biased than the square sampling of image texture. This is formalized in information theory where hexagonal structures are shown to be optimal quantizers over 2D regular lattices [START_REF] Newman | The hexagon theorem[END_REF]. Second, there is no discontinuity in the appearance model, nor any distortion. This makes filtering more accurate and much easier to perform in practice, as illustrated in Figure 5. Third, sampling frequency can be chosen locally, at the face level, whereas it is fixed globally with image textures (local sampling frequency depends on the vertical and horizontal texture resolution and the mapping function). This is particularly suitable for multi-view 3D modeling, as detailed in the next section. Fourth, local editing or resolution changes do not require a complete resampling or recharting of the mapping.

Figure 5: Illustration of a filtering operation with a texture from real images. Three close-ups of a mesh, rendered with: Left: the unfiltered image texture, middle: the same image texture with a simple sharpening filter applied, right: the mesh texture with a sharpening operation applied directly on the mesh. Artifacts appear despite the texture being dilated many times in the image texture. They are due to the seams lying too close to peaks in the color gradient.

Storage and Compression.

Similar to [24], we decouple geometry and color information, to make it suitable for standard graphics pipeline. Each vertex, edge and face stores an index to a single mesh global color array. No compression scheme is yet available for mesh textures and we introduce a novel approach for that purpose in the next section 4.2.

Mesh Textures from Real Images

As previously discussed, mesh textures are convenient when modeling shapes using real images. We detail in this section the related aspects in that case, namely the strategies for the appearance sampling and the compression.

Sampling Strategy

Given the pixel information in the observed images we define an adaptive sampling strategy that optimally exploits the appearance data. To this aim, we first choose a resolution level per triangle before computing the mesh texture. Then, in a post-processing step, we downsample triangles when they can still be accurately interpolated by the next lower resolution level, this up to a given error threshold. The first step depends directly on the pixel density in the input images whereas the second step depends on the amount of information present in the pixels. Figure 6 illustrates the benefits of this sampling strategy: The uniform sampling of texels with the image texture implies that unseen or uniform areas are oversampled, while others could benefit from a denser sampling, e.g. the button on the coat. In contrast, the mesh texture representation samples fewer points in the hidden parts of the mesh but has a denser sampling in this highly textured area.

Adaptive Sampling

Our sampling strategy adapts the resolution R i of the triangle T i to the number of pixels that observe this triangle. To this aim, we assume that, for any subpart S ⊂ M of the mesh, a function D gives a measure of the available input pixels relevant to S, dependent on the coloring strategy chosen. For example, a super-resolution strategy can define D as the total number of pixels, in all selected views, that project onto S. A naive blending approach will define D as the maximum number of pixels, among all views, that project onto S. In this paper, we implement a strategy that selects a single best camera per triangle (see 5.1). Thus, D measures the number of pixels in this selected view. Now given the number D(T i ) of pixels in triangle T i we need to find the appropriate appearance resolution R i for that triangle. For such a resolution, the number of appearance samples in the triangle is approximately:

(R i -2)(R i -1) 2 + 3 2 (R i -1) + 3 6 = 1 2 R 2 i , (4)
counting samples within the triangle as full, samples over edges as shared with 2 triangles and samples on vertices as shared, on average, among 6 triangles. Given this number, for each triangle T i , we choose the smallest possible resolution R i that is a power of two and such that:

1 2 λR 2 i < D(T i ) , (5)
where the parameter λ can be chosen depending on the ex-pected rendering quality. Typically, λ = 2 in our experiments.

Downsampling

Real world objects often present large regions with more or less uniform appearances. In such a region, irrespective of the resolution R i , the appearance samples of a triangle will carry redundant color information. In order to account for that, we downsample triangles that fall into this category. More precisely, for each triangle T i , we consider all sample points associated to the current resolution level and compare their color values with the ones they would get by just interpolating the next lower-resolution level. Looking back at equations ( 2) and ( 3), these are samples with barycentric coordinates m odd, or n odd, or both. Denoting o = R -m -n, we compute the average difference as:

E = 8 3n(n + 1) × (Σ m + Σ n + Σ o ) , (6)
with:

Σ m = m even, n,o odd d(s (m,n,o) , s (m,n-1,o+1) + s (m,n+1,o-1) 2 ) , (7)

Mesh Texture Compression

One of the important features with image textures is the ability to apply efficient compression schemes and to significantly reduce the amount of data while keeping good visual quality. In order to make mesh textures an efficient appearance modeling tool as well, we also give them this ability.

Strategy

Image compression is a widely studied problem for which optimal solutions have been proposed. Since mesh textures store color information in 2D (triangular) lattices, we believe that the existing schemes already provide efficient solutions with, furthermore, standard components. Most state-of-the-art lossy image compression and texture compression methods, including JPEG, ASTC or ETC, follow a block-based approach, first introduced by Delp and Mitchell with Block Truncating Coding in 1979 [5], in which they decompose the image into rectangular blocks and encode each block independently. We draw inspiration from JPEG [20] and adapt it to mesh textures. We first give an outline of our method below. The full process is then detailed in the next subsection (4.2.2).

Vector space -The triangles of a mesh texture are not 2D regular grids but they constitute anyway natural 'block' candidates. However, block sizes are predetermined and cannot be chosen arbitrarily: we need to be able to process blocks of varying size. Since we only allow resolution levels that are powers of two, the number of potential block sizes is limited (typically 4 or 5). We process each block size (i.e. each vector space) independently.

Space transform -In our case, blocks are triangular rather than rectangular and we cannot apply DCT directly. While solutions have been proposed to adapt DCT to triangular shapes (e.g. [6]), we opted for a different approach based on PCA decompositions, this after having tested both methods. While this transformation has little physical meaning compared to DCT, it is an effective tool for identifying the few components that encompass most of the variations in the data. Besides, we believe the main components also tend to cover low-frequency variations that are most relevant to the human visual system. (See figure 7).

Discarding irrelevant components -Similar to JPEG, we downsample chroma components. Coefficients in the PCA basis are quantized based on a specific quantization matrix.

Entropy coding -Huffman coding and run-length encoding: We choose to compute Huffman tables for each case (i.e. mesh), rather than use predefined ones, as in standard JPEG. We compute separate Huffman tables for each resolution level. The PCA decomposition gives us a natural writing order for coefficients (as opposed to the zigzagging needed for JPEG).

Encoding and Implementation.

We present here the full encoding pipeline in more details. For our compression scheme, edge colors are duplicated. While this introduces redundancy in the data, it allows us to compress edges' color information within the 2D signal of the triangles. More importantly, it completely eliminates the need for indexing between geometry and color data. Vertex colors are stored separately and left uncompressed. The encoding process is as follows:

• Triangles are reordered per resolution level R, from highest to lowest. Each set of triangles with a given resolution is processed independently.

• RGB data is transformed into Y C b C r . Chroma components are downsampled to the next lower resolution. This means eliminating between 1/2 (if R = 2) and 3/4 (asymptotic behaviour when R increases) of chroma samples.

• A PCA decomposition is computed for each value of R.

• The coordinates of faces in this new space are quantized. We choose quantization vectors of the form F loor(1 + ai + bi 2 times a normalization factor, where i is the PCA component index. i ranges from 0 for the main component, to (D T -1), and a and b are chosen empirically (in our experiments, a = 1, and b ∈ [0.01, 0.1] for varying bitrates).

• The mean vector and eigen vectors of the PCA decomposition are also quantized (e.g. on 12 bits per coefficient).

• Barring some implementation details, we then follow the standard JPEG pipeline, except that coefficients are not subdivided into several channels. Resolution changes are written using a specific JPEG marker. They are followed by the quantization table and the list of PCA eigenvectors (and mean), signaled by their own marker as well. Then,we write the Huffman tables used for this part of the data, and finally, the entropy-coded data.

The geometry data is written in another file. It includes a list of vertices (3D points) with one color per vertex, and a list of faces, i.e. three vertex indices. Faces are written following the order used in the compressed color file.

The color data of each face takes up a varying number of bytes in the encoded signal, which cannot be predicted. Thus, for rendering, data must be decoded before being loaded on memory (which is also the case with JPEG textures, but not with ETC or ASTC, for example ( [START_REF] Ström | i PACKMAN: high-quality, low-complexity texture compression for mobile phones[END_REF]14]).

Evaluation

In order to demonstrate the benefit of mesh texture in the context of multi-view shape reconstruction, we report in this section on comparisons between image and mesh textures given reconstructed shapes and their observed images. The criteria we consider are visual quality or fidelity to the original images and size which, we believe, are the most critical properties of the appearance in our context.

We first detail the appearance function that associates a color to any point on the shape surface (section 5.1).. This function is used with both image and mesh textures in the comparisons. We then render views of our models for comparisons, both in terms of sampling (section 5.2) and compression ratio (section 5.3) and with respect to the criteria mentioned. Additional experiments with more data are included in the supplementary material.

Appearance Function

Given a mesh and its associated observed images from different viewpoints, we compute a continuous function that gives a color value for each point on the surface mesh. Various strategies can be considered for that purpose, from simple blending to super-resolution approaches. Our objective is primarily to compare appearance representations and not to evaluate appearance function, we opt therefore for an effective best view strategy that select for each mesh face the best view available. The key aspects are as follows:

Computing visibility: Each input image I i is upsampled, and visibility for every (face,view) pair is computed at the sub-pixel level, which gives more precision on the available appearance information.

Assigning Views: In real acquisition setups, the camera properties are inaccurate and the input views are consequently misaligned and inconsistent. Averaging pixel information over different images is therefore likely to severely blur the appearance. Following many works in that respect, e.g. [11, 7, 22], we assign a single view to each face of the shape model. In practice, we select the most informative view with the largest number of visible subpixels.

Blending views: The view assignment can be seen as weighing the input views by a function ω 0 that depends on the position on the surface and the view index. Choosing a view j for face T means

∀i ∈ {1, ..., I}, ∀p ∈]0, 1[ 2 , ω 0 (T , p, i) = 1 if i = j 0 otherwise . (8)
To avoid visible seams between triangles, we compute an estimate of ω 0 for each vertex by averaging its value on adjacent faces. Weights are then interpolated smoothly between vertices, i.e. for the image i, a face t = {v 1 , v 2 , v 3 }, a location p = (a, b, c) within t, we use the function ω defined by:

ω(x(t, p), i) = a ω 0 (v 1 , i)+b ω 0 (v 2 , i)+c ω 0 (v 3 , i) . (9)
Finally, the color C of each surface point is defined by where {I 1 . . . I I } are the input images with their respective 3D-2D projection operators {π i }.

∀x ∈ M, C(x) = I i=1 I i π i (x) × ω(x, i) , (10)

Sampling

We evaluate our method on real data, captured with multi-camera platforms. The input data and the reconstructed mesh are provided by the authors of [START_REF] Leroy | Shape reconstruction using volume sweeping and learned photoconsistency[END_REF]. It provides 64 different input views of the same scene, and includes an image texture, computed with a method based on conformal maps. We use a leave-one-out evaluation strategy, removing a given input view, and comparing it against the projection of appearances computed using the remaining views. We sample the function C defined in equation (10) with varying parameters. Charts in image textures are dilated to prevent artifacts and for a fair comparison with mesh textures. We use two different metrics for such comparison: the Multi-scale Structural Similarity (MS-SSIM) from [START_REF] Wang | Multiscale structural similarity for image quality assessment[END_REF], and the more recent Learned Perceptual Image Patch Similarity (LPIPS) from [START_REF] Zhang | The unreasonable effectiveness of deep features as a perceptual metric[END_REF]. We compare images within the mesh silhouette only in the images. Figure 8 shows some numerical results. MS-SSIM measures a similarity, thus higher scores are better, contrary to LPIPS which measures the perceptual difference. In this case, MS-SSIM seems ill-suited for a meaningful analysis, given the very small range of variation. At high resolution, sampling is not a limiting factor, and the error is mostly due to the imperfect geometry and color function C. As we progressively decrease the number of samples, LPIPS shows our representation tends to retain more visual information.

Compression Evaluation

To evaluate our compression method, we pick a highresolution texture image and texture mesh of the same mesh, with a similar number of samples and a similar score on both metrics. We compress them with varying quantization matrices (or compression ratio for JPEG), and compare the results, using the same process as in the previous section, except that this time, we use a single model that we backproject on different input views. Figure 9 demonstrates that our representation outperforms the image texture, especially at low bitrates. Figure 2 displays the artifacts inher- ent to both methods. If we zoom in, we can already notice small discontinuities with the image textures, because of seams in the texture in that case. As we compress the texture with increasingly high compression ratios, blocky artifacts start to appear along the seams, and finally, on the whole mesh. By comparison our mesh texture method yields triangular blocky artifacts that seem less perceptible, probably because they follow an irregular pattern. Besides, our PCA decomposition basis is computed specifically for the mesh, contrary to the more general DCT decomposition.

Conclusion

We have studied the benefits of representing color information directly on 3D meshes in the context of multiview image-based appearance modeling. More specifically, we have compared the use of our proposed mesh texture pipeline, with the more widespread use of texture maps. First, taking advantage of its sampling flexibility over the surface triangles, we were able to formulate a smart locallyadaptive sampling strategy, which we show to be more efficient than uniform sampling in the image-based appearance modeling context. Second, we introduce a novel compression strategy dedicated to our mesh-embedded color information, and show that it is effective enough to outperform classic texture storage in terms of the perceptual fidelity vs storage capacity tradeoff, especially at low bit rates. The results obtained validate this scheme as a practical solution, even when dealing with large amounts of data. For future work, it would be interesting to implement compression strategies with a fixed bit rate, so that meshes could be rendered directly from their compressed mesh texture. Additionally, one could try to build temporally consistent mesh textures, to represent the appearance of a moving object across several frames in a compact way. Finally, we would like to combine this framework with a state-of-theart superresolution algorithm, to see if additional quality performance can be achieved with this representation.

Figure 2 :

 2 Figure 2: Visualization and surface close-up of our appearance compression scheme: Top: mesh texture. Bottom: image texture. Left to right: high-resolution texture, texture compressed to decreasing file size.

Figure 3 :

 3 Figure 3: Close-up of the image texture of a mesh, generated with the same texture atlas but different resolutions. (Left: 80 × 50, Right: 320 × 200).

Figure 4 :

 4 Figure 4: Appearance samples on a mesh texture model with triangles exhibiting varying appearance resolutions: from left to right, R = 8, 4, 2, 1.

Figure 6 :

 6 Figure 6: Right: Two close-up renderings of the same colored mesh, with an image texture (top) and a mesh texture (bottom), using approximately the same number of samples, rendered with closest point interpolation. (Close-ups of the model on the left)

  where s (m,n,o) is the color of the sample with barycentric coordinates (m, n, o) and d is the L 2 distance in corresponding color space. If E is less than a threshold T ds , we downsample the triangle appearance and verify again at the next lower resolution. In our experiments, T ds ∈ [0, 80] for varying visual quality (with color values in [0, 255]).

Figure 7 :

 7 Figure 7: Example of the PCA decomposition for triangles of resolution 8 for one of our test meshes. Each triangle represents one component, from most to least relevant (i.e. biggest to smallest eigenvalue) in reading order. The first few components (which are most preversed by the quantization step) encode low frequency information.

Figure 8 :

 8 Figure 8: Result on the test mesh: score averaged over all viewpoints. (red) Image textures, (blue) Mesh textures.

Figure 9 :

 9 Figure 9: Result on the test mesh: score averaged on different viewpoints. Image textures shown in red, mesh textures in blue. The original image texture takes 5.6 MB (3072 × 3072) for 2940k texels. The original mesh texture takes 5.5 MB, for 3007k color samples.


Carole Duffet
email: duffet@math.univ-montp2.fr

Delphine Sinoquet

Quantifying uncertainties on the solution model of seismic tomography

Reflection tomography allows the determination of a propagation velocity model that fits the traveltime data associated with reflections of seismic waves in the subsurface. A least-squares formulation is used to compare the observed traveltimes and the traveltimes computed by the forward operator based on a ray tracing. The solution of this inverse problem is only one among many possible models. A linearized a posteriori analysis is then crucial to quantify the range of admissible models we can obtain from these data and the a priori information. The contribution of this paper is to propose a formalism which allows us to compute uncertainties on relevant geological quantities for a reduced computational time. Nevertheless, this approach is only valid in the vicinity of the solution model (linearized framework), complex cases may thus require a nonlinear approach. Application on a 2D real data set illustrates the linearized approach to quantify uncertainties on the solution of seismic tomography. Finally, the limitations of this approach are discussed.

Introduction

Geophysical methods for imaging complex geological subsurfaces in petroleum exploration requires the determination of an accurate propagation velocity model. Seismic reflection tomography (see [START_REF] Bishop | Tomographic determination of velocity and depth in laterally varying media[END_REF]) turns out to be an efficient method for that: this method allows us to determine a seismic velocity distribution from traveltime data associated with the seismic waves reflecting on geological surfaces. This inverse problem is formulated as a leastsquares minimization problem which consists of the minimization of the mismatch between the observed traveltimes and the traveltimes computed by the forward problem (solved by a ray tracing method). The first section of this paper presents an overview of the reflection tomography method we have developed (see [START_REF] Jurado | 3D reflection tomography designed for complex structures Expanded Abstracts[END_REF][START_REF] Jurado | Fast 3D two-point ray tracing for traveltime tomography[END_REF][16][START_REF] Stopin | Reflection tomography of compressional and shear modes for determination of anisotropic velocity models[END_REF] for more details).

In principle, the estimation of the uncertainties begins with an analysis of the errors on the solution model by applying Monte Carlo techniques among others (for examples see [START_REF] Mosegaard | Monte Carlo sampling of solutions to inverse problems[END_REF][START_REF] Malinverno | A Monte Carlo method to quantify uncertainty in the inversion of zero-offset[END_REF]). These methods have a prohibitively expensive computational cost. Our goal here is to develop methods for the a posteriori analysis of the solution which will be easier to apply and nevertheless provide useful results. Classically, the analysis of the solution consists only of checking the misfits between observed traveltimes and calculated traveltimes by computing characteristic values of the misfit distribution and by studying the spatial distribution of the misfits. This analysis is not sufficient: even if these quality control criteria are matched, the determined model is only one of many possible models that match the data. An uncertainty analysis should be performed to quantify the range of admissible models we can obtain from these data and the a priori information (see [START_REF] Duffet | Quantifying geological uncertainties in velocity model building Expanded Abstracts[END_REF]).

The different methods we are interested in here allow us to access the uncertainties on the solution model thanks to the analysis of the a posteriori covariance matrix obtained in the linearized framework (see [START_REF] Gouveia | Resolution of seismic waveform inversion: Bayes versus Occam Inverse[END_REF][START_REF] Tarantola | Inverse Problem Theory[END_REF] among others). The computation of this matrix is generally expensive for 3D problems; among others, [START_REF] Berryman | Analysis of approximate inverses in tomography: II[END_REF] proposed a method to compute approximations of the covariance and resolution matrices from the different iterates when the Gauss-Newton linear system is solved by an iterative method such as the gradient conjugate algorithm. But the physical interpretation of the terms of this huge matrix may be cumbersome when we want to go further than the classical but restrictive analysis of the diagonal terms. In this paper, we propose another approach which simulates admissible models from an a posteriori probability density function that directly provides physical models. This method requires a Cholesky decomposition of the Gauss-Newton matrix which may be unfeasible for 3D realistic models. To overcome this difficulty, we propose a formalism based on linear combinations of model parameters (macro-parameters) that allows us to compute uncertainties on relevant geological quantities (the average thickness of a layer for example) for a reduced computational time (the a posteriori covariance matrix is reduced to the macro-parameter space). The application of these methods on two real data sets shows their effectiveness and their limitations are discussed, particularly the limitations due to the linearization: a nonlinear approach is proposed in order to delimit the space of admissible solutions.

In this paper, after the description of the chosen formulation of reflection tomography in section 2, we present the two proposed methods in section 4 to perform a linearized uncertainty analysis on the solution model of a 2D real data set described in section 3. In section 5, we explain the motivations to carry out a nonlinear approach on this 2D complex case and present some results. Finally, section 6 is devoted to concluding the paper.

The reflection tomography problem

Let us first present the reflection tomography method which was developed in the KIM consortium [16]. The chosen model representation is a blocky velocity model where the velocity distribution is described by slowly varying layer velocities (also called velocity blocks) delimited by interfaces (see figure 1 for an example). The model is thus composed of two kinds of parameters: parameters describing the velocity variations within the layers and parameters describing the geometry of the interfaces delimiting the velocity blocks 3 . Moreover, anisotropy of the velocity may be modelled by two parameters η and δ (see [START_REF] Stopin | Reflection tomography of compressional and shear modes for determination of anisotropic velocity models[END_REF] x(km)

x for more details). This blocky representation allows the representation of complex models (e.g., presenting rapid lateral velocity variations and complex geological structures) with a reasonable number of unknowns. Strong a priori information is introduced, the number of layers being fixed: this reduces the underdetermination of the inverse problem.

The forward problem of reflection tomography consists of computing traveltimes of reflected seismic waves, given a model, an acquisition survey (locations of the sources and the receivers) and signatures (generally reflectors where the waves reflect). It is solved by a ray tracing method which is a high frequency approximation of the wave equation (see [START_REF] Cervený | Raytracing algorithms in three-dimensional laterally varying layered structures Seismic Tomography[END_REF][START_REF] Jurado | Fast 3D two-point ray tracing for traveltime tomography[END_REF] among others). We denote by T (m) the forward modelling operator, that gives, for a specified reflector, the traveltimes associated with all source-receiver pairs. Reflection traveltime tomography is the corresponding inverse problem: its purpose is to adjust m such that T (m) best matches a vector of traveltimes T obs picked on seismic data. A natural formulation of this problem is the least-squares formulation

T (m) -T obs 2 C -1 D + m -m prior 2 C -1 M , (1)
where

• -A denotes the norm defined by A matrix: x A = x T Ax.

• C D is the a priori covariance operator in the data space that allows us to describe errors on the data, • m prior is an a priori model (coming, for instance, from geological knowledge and/or additional information resulting from well measurements), • C M is the a priori covariance matrix in the model space that describes errors on the a priori model.

In practice, building an a priori model is difficult and the a priori term, which is necessary to well pose the inverse problem, is often replaced by regularization terms made up of second derivatives of the model (see for instance [START_REF] Lailly | Ill-posed and well-posed formulations of the reflection traveltime tomography problem[END_REF]) which leads to

C -1 M = Q reg and m prior = 0 T (m) -T obs 2 C -1 D + 2 D 2 m 2 Id (2)
where a continuation technique (see [START_REF] Bube | A continuation approach to regularization for traveltime tomography[END_REF]) is applied for tuning the regularization weight : the weight is decreased until the model matches the data with the expected accuracy. The obtained model is then the smoothest model which fits the data. This methodology allows us to stabilize the inversion. This large size (up to ∼100 000-500 000 traveltime data and ∼5000-10 000 unknowns) nonlinear least-squares problem (the objective function is not quadratic since the forward modelling operator T (m) is nonlinear) is solved classically by a Gauss-Newton method based on successive linearizations of the forward operator which needs the computation of the Jacobian matrix of T J (m) = ∂T ∂m (m) . 4 The resulting approximation of the Hessian matrix of the cost function is noted H (m) = J (m) t C -1 D J (m). Then, the quadratic approximation of the objective function ( 2)

Cn (δm) = J n δm -δT obs n 2 C -1 d + 2 D 2 (m n + δm) 2 (3)
is minimized at each Gauss-Newton iteration n by a preconditioned conjugate gradient (m n is the current model, δm is the unknown model perturbation, J n = ∂T ∂m (m n ) is the Jacobian matrix evaluated at m n and δT obs n = T obs -T cal (m n )) (see [START_REF] Delbos | Constrained optimization in seismic reflection tomography: a Gauss-Newton augmented Lagrangian approach[END_REF]).

Presentation of the application

Let us consider a 2D real data set already studied by [START_REF] Broto | Anisotropic traveltime tomography for depth consistent imaging of PP and PS data The[END_REF] with PP and PS 5 data to illustrate all the following. 45 338 traveltime data were interpreted and an uncertainty of 5 ms (respectively 8 ms) is associated with PP data (respectively PS data). A layer-stripping approach (separate inversion of each velocity layer from the shallowest layer to the deepest one) provided the velocity model of figure 1.

The model is described by four interfaces, corresponding to the interpreted events (H 1, H 3, H 4 and H 5) which define three layers with only lateral velocity variations for the first two uppermost layers (v(x), there is no vertical variations) and the deepest layer stretching from H 3 to H 5 with a 2D velocity VP5(x, z) and VS5(x, z) (figure 1). VP5 and VS5 define respectively the vertical velocity propagation of the P-waves and the velocity propagation of the S-waves. This model is composed of 4588 parameters, 592 for the interfaces and 3936 for the velocities. In the deepest layer, we assume a vertical transverse isotropic velocity field characterized by VP5 and VS5 and the two parameters η and δ measuring respectively the measure of the velocity an-ellipticity and the near vertical velocity propagation of the P-waves 6 . The values of η and δ in the model correspond to a strongly anisotropic medium. For the first two layers, the inversion results are very satisfactory, with a traveltime misfit RMS (root mean square) of 3.5 ms for the PP-data and of 7 ms for the PS-data. For the last layer, the traveltime misfit RMS is of 3.7 ms and of 5.6 ms for, respectively, H 4 and H 5 PP-data and of 8.5 ms for H 5 PS-data, results which are consistent with the data uncertainties (figure 2).

As already shown by [START_REF] Stopin | Reflection tomography of compressional and shear modes for determination of anisotropic velocity models[END_REF], it turns out that the anisotropy parameter δ is strongly undetermined from seismic data. The value of δ parameter was obtained in [START_REF] Broto | Anisotropic traveltime tomography for depth consistent imaging of PP and PS data The[END_REF] by a trial and error approach in order to match approximately the depth of H 5 horizon given at well location. We propose to carry on an a posteriori uncertainty analysis in order to quantify the uncertainties on this solution model: we focus on the anisotropic layer delimited by H 3 and H 5 (velocities VP5 and VS5).

Quantifying uncertainties on the solution model: a linearized approach

A classical approach to quantify uncertainties consists of the analysis of the Hessian matrix (or its inverse: the a posteriori covariance matrix) associated with the linearized problem (3) around the solution m ∞ . This approach is only valid in the vicinity of the solution model, the size of the vicinity depending on the nonlinearity of the forward map.

The bi-linear form associated with the Hessian matrix measures the influence of a model perturbation δm on the quadratic cost function defined around the solution m ∞ :

C(δm) -C(0) = 1 2 (J ∞ δm) T C -1 D (J ∞ δm) + 1 2 δm T C -1 M δm -δm T J T ∞ C -1 D δT obs (m ∞ ) -C -1 M m ∞ = 1 2 δm T J T ∞ C -1 D J ∞ + C -1 M δm = 1 2 δm T H (m ∞ )δm. (4) with • J ∞ = J (m ∞ ), • C(0) = C(m ∞ ), and • J T ∞ C -1 D δT obs (m ∞ ) -C -1 M m ∞ = g(m ∞ ) = 0, gradient of C(m) that vanishes at the solution m ∞ .
The a posteriori covariance matrix is defined by the inverse of H (m ∞ ) (see [START_REF] Tarantola | Inverse Problem Theory[END_REF])

C M = J T ∞ C -1 D J ∞ + C -1 M -1 . ( 5
)
The space of admissible models can be characterized by the contour lines

(m -m ∞ ) T C -1 M (m -m ∞ ) = constant, (6)
which are ellipsoids of centre m ∞ and correspond also to contour lines of the a posteriori Gaussian probability density function

exp -1 2 δm T C -1 M δm . ( 7
)
The diagonal terms of C M are the uncertainties on the parameters describing the model and the off-diagonal terms are the correlations between these uncertainties. For instance, the probability that the true model parameter

p i verifies -2(C M ) i,i p i -p ∞ 2(C M ) i,i
, independently of the values of the other model parameters, is about 95%. To take into account the correlations between the parameters, we should study the 95% confidence ellipsoid. The axes of the ellipsoids (6) are defined by the eigenvectors of C M , the square root of the eigenvalues giving the uncertainties on the associated eigenvector.

Remark 1. In practice, the diagonal terms of C M provide the uncertainties on the B-spline parameters, unknowns of the discretized inverse problem: they are not physical quantities. We would rather compute uncertainties on physical quantities, for instance, the evaluation of the B-spline functions in the physical domain.

Remark 2. Following [START_REF] Lailly | What information on the earth model do reflection travel times hold?[END_REF], the eigenvector decomposition of the a posteriori covariance matrix gives access to the worst/best determined model directions and the associated uncertainties. The best (respectively worst) determined model direction corresponds to the smallest (respectively highest) eigenvalue of the a posteriori covariance matrix, i.e. highest (respectively smallest) perturbation of the quadratic cost function. The uncertainties associated with these eigenvectors, namely the square root of the eigenvalues, are meaningless, the eigenvectors being composed of mixed velocity and interface parameters, combinations which are usually difficult to link to physical quantities.

From those remarks, we propose in the next sections two methods to quantify geological uncertainties on the solution model which avoid the expensive computation and the cumbersome analysis of the generally huge a posteriori covariance matrix C M .

Simulations of admissible models

The first proposed method to quantify uncertainties is the simulation of admissible models from the a posteriori probability density function (pdf) [START_REF] Delbos | Constrained optimization in seismic reflection tomography: a Gauss-Newton augmented Lagrangian approach[END_REF]. The simulations provide directly interpretable results, i.e. physical models. The method (see for instance [START_REF] Parker | Geophysical Inverse Theory[END_REF]) consists in random simulations of model perturbations following the pdf [START_REF] Delbos | Constrained optimization in seismic reflection tomography: a Gauss-Newton augmented Lagrangian approach[END_REF].

First, we apply the variable transformation:

δm = Uδm, ( 8
)
where U is the lower triangular matrix obtained by a Cholesky decomposition of the inverse of the a posteriori covariance matrix C -1 M = H = U T U (symmetric positive definite matrix). The method consists then in simulations of the uncorrelated Gaussian pdf with unit variance exp -1 2 δm T δm .

Simulations δm are transformed into correlated Gaussian simulations of vector δm by the inverse transformation of (8). Figures 34567show 100 simulated models obtained by this method from the solution model of figure 1: interfaces H 4 and H 5, variations of velocities VP5 and VS5 along x and z directions and histograms of the anisotropy parameters η and δ. From these simulations, we observe that the highest uncertainties on the lateral velocity variations are located at the boundaries of the model, areas that are not well illuminated by the rays. Indeed, for a slice at constant z = 2.7 km, we observe uncertainties of ≈690 m s -1 for VP5 and of ≈260 m s -1 for VS5 at the boundaries of the model and uncertainties of ≈450 m s -1 for VP5 and of ≈140 m s -1 for VS5 in the illuminated parts of the model. Concerning the interface depths, we observe uncertainties of ≈180 m for reflector H 4 and of ≈220 m for reflector H 5 at the boundaries and of ≈90 m for reflector H 4 and of ≈100 m for reflector H 5 elsewhere. For the anisotropy parameters η and δ, we notice uncertainties around 0.3% on η and 2% on δ.

This method is quite attractive for the straightforward interpretation of the results despite its cost (cost of the Cholesky decomposition): it provides physical models leading to small perturbations of the quadratic cost functions (3).

Uncertainties on geological macro-parameters

We propose here a method that allows us to deal with large size models at a reasonable cost and provides uncertainties on chosen physical quantities. The proposed approach consists in building macro-parameters (MP) with a geophysical interest. These macro-parameters are linear combinations of the inverted parameters such as the mean of the velocity variations in a zone, the slope of an interface, the average thickness of a layer, etc. [START_REF] Grenié | Tomographie de temps de trajet adaptée au problème de l'imagerie 3D en exploration pétrolière: théorie, mise en oeuvre logicielle et applications[END_REF] has introduced the notion of macro-parameter (his main motivation being to avoid numerical problems in the inversion of the complete Hessian). We propose here a generalization of his work (general definition of macro-parameters) which allows the computation of uncertainties for huge 3D problems.

We define a macro-parameter as

P = Bp, ( 10
)
where p is the n p model parameters vector (B-spline parameters in our case), P is the n MP macro-parameters vector and B is the condensation matrix. We compute the a posteriori covariance matrix C M in the macro-parameter space:

C M = BC M B T . ( 11
)
Note that C M (an n MP × n MP matrix) is small compared to C M (an n p × n p matrix) since n MP n p . To obtain this matrix, we do not need to compute the whole inverse of the Hessian, C M = H -1 , in the parameter space.

Indeed, to obtain C M = H -1 = BH -1 B T in the MP space, we just need H -1 B T j , where B T j are the different columns of B T . We thus solve n MP linear systems, H -1 j = BH -1 B T j , that are similar to the linearized problem we solve at each Gauss-Newton iteration. Thus, the computational cost for one MP is comparable to one iteration of the inversion process.

We applied this method on the model of figure 1. We choose simple MP: the mean of the velocity variations and the mean of the interface depth in the illuminated part of the layer (region covered by rays) 7 . The results are listed in table 1. The uncertainty on anisotropy parameter δ is high: ±74% (if we consider twice the relative standard deviation) and we observe high correlations between η and δ (0.93) and also correlations between anisotropy parameters and H 5 depth (-0.33 and -0.36). We notice also the bad determination of the velocities: ±22% (relative value). All these results are consistent with the simulation results (section 4.1).

This method with its general formalism allows us to compute uncertainties on relevant geological quantities with a reasonable computation cost and to highlight easily strong coupling between the chosen MP.

Discussions on the two methods

As already mentioned, the two methods described in sections 4.1 and 4.2 rely on the quadratic approximation (3) of the nonlinear cost function. In table 2 we have listed the RMS of traveltime misfits for 20 simulated models (computed by the ray tracing forward operator). Note that the RMS traveltime misfit for the solution model is 6.2 ms. Thus, there is only one model which is not acceptable (11.6 ms for the model 15). The others are admissible models with a RMS traveltime misfit bounded by 7 ms. The simulation method has then provided model perturbations that correspond to small perturbations of the quadratic cost function but also to small perturbations of the nonlinear cost function (RMS of traveltime misfits remains small). It shows for this example, the good agreement between the quadratic cost function and the nonlinear one around the solution model. Figure 8 illustrates isovalues of the probability density function [START_REF] Delbos | Constrained optimization in seismic reflection tomography: a Gauss-Newton augmented Lagrangian approach[END_REF] in the MP space, i.e. the ellipsoids defined by δm T MP C -1 M δm MP . We visualize the 2D-marginal probability density functions for several relevant MP couples, the 68% confidence ellipsoid and the 95% confidence ellipsoid for η and δ, VP5 and VS5, δ and H 5. This representation allows us to underline the correlations between the anisotropy parameters η and δ which are quite strong in this example. The 68% and the 95% confidence ellipsoids are elongated (the left part of figure 8). On the opposite, the P and S velocities are uncorrelated: we obtain almost a circle (the middle part of figure 8). We can also observe correlations between the anisotropy parameter δ and the mean of the last reflector depth H 5 (the right part of figure 8). In this figure, we have also superimposed the anisotropy parameters η and δ and the MP VP5 and VS5, δ and MP H 5 built from the 100 simulated models. This figure allows us to point out the equivalence of the simulation method and the macro-parameter approach.

Remark that the macro-parameter approach provides interesting information on the correlations of macro-parameter uncertainties and has a lower computation cost than simulations of the model parameters using the complete a posteriori covariance matrix. We could also perform Gaussian simulations of the macro-parameters using the reduced a posteriori covariance matrix.

Nonlinear a posteriori analysis

We have shown in the previous sections that a linearized uncertainty analysis allows us to delimit the range of possible solution models that fit, with the expected accuracy, the data and the a priori information. Nevertheless, we should keep in mind that this approach is only valid in the vicinity of the solution model (linearized framework) and complex cases may require a nonlinear approach.

Model 15 of table 2 indicates that some of the simulated models may produce unacceptable traveltimes misfits (RMS ≈ 11 ms and MAX = 25 ms for PP data and =41 ms for PS data). It shows the limitation of the quadratic approximation of the nonlinear cost function and thus the limitation of the linearized methods of uncertainty quantification. Moreover, [START_REF] Delbos | Constrained optimization in seismic reflection tomography: a Gauss-Newton augmented Lagrangian approach[END_REF] has shown that a tomographic inversion with constraints on the location of the interfaces at well locations provided a very different model (model 2 shown in figure 9) from the solution model of figure 1. The two models verify the traveltime data with the expected accuracy. This second model does not belong to the range of possible models detected by the two methods we proposed: indeed, δ is equal to 15.94% whereas the expected range of values was between -8% and -1% (see for instance figure 7). It shows again the limitations of the linearized approach. To perform a nonlinear analysis we have chosen an experimental approach which consists of performing several constrained inversions allowing us to test different geological scenarii to try to delimit the space of admissible solutions. A sequential quadratic programming (SQP) approach based on an augmented Lagrangian method (for details see [START_REF] Delbos | Constrained optimization in seismic reflection tomography: a Gauss-Newton augmented Lagrangian approach[END_REF]) is used to solve this nonlinear constrained optimization problem (computational times listed in table 3).

For instance, we could test different hypotheses on the values of the anisotropy parameters for which a strong uncertainty has been detected by the linearized approach. If we introduce a constraint on the anisotropy parameter δ, such as δ > 0, we find model 3 displayed in figure 10: this result expands the range of admissible models (variability of δ detected by the MP approach was -4.43 ± 3.2%).

In the same way, we can test the variability of the anisotropy parameter η. We note (see table 4) that the value of the anisotropy parameter η of model 2 (figure 9) is 8.27%, but all the simulated values of this parameter are around 6% and the variability of η provided by the MP is ±0.44%. Nevertheless, tomographic inversion under the constraint η 4% or η 8.5% does not provide a model that satisfies both these constraints and the data with the expected accuracy.

We have shown that, by constrained tomography, we can test different geological hypotheses and delimit the space of admissible models and go further than the linearized approach which explores only the vicinity of the solution model which may be one local minimum around others.

Conclusions-discussions

Reflection tomography provides the velocity model that best fits the traveltime data. However, this solution is only one among many admissible models. An a posteriori uncertainty analysis is crucial to delimit the range of possible solution models that would fit the data and the a priori information with the expected accuracy. In this paper, we describe two methods, the simulations and the macro-parameter approach, to perform a linearized a posteriori analysis, approach valid only in the vicinity of the solution model. The simulation method based on the analysis of the a posteriori covariance matrix is quite attractive for the straightforward interpretation of the results despite its cost. We propose a general formalism to reduce the a posteriori analysis to geological quantities of the model and thus reduce the cost of the uncertainty estimation. These methods give interesting results on a 2D real data set. These methods being valid only in the vicinity of the solution we propose to complete the linearized study by a nonlinear approach based on constrained optimization.

Figure 1 .

 1 Figure 1. Solution velocity model obtained by tomography. Left: P velocity. Right: S velocity. The RMS value of the traveltime misfits is 6.2 ms. The anisotropy parameters values are η = 6.29% and δ = -4.43%.

Figure 2 .

 2 Figure 2. Interfaces of the solution velocity model obtained by tomography and their associated traveltime misfits: impact points of rays reflected on H 4 and H 5 (Left: PP rays. Right: PS rays), scales represent traveltime misfits. Note that the illuminated zone extends from 14 km to 20.5 km.

Figure 3 .Figure 4 .Figure 5 .

 345 Figure 3. Reflector H 4 extracted from the 100 simulated velocity models. The bold line represents the solution interface.

Figure 6 .Figure 7 .

 67 Figure 6. S-velocity variations along the x axis at z = 2.7 km (left) and along the z axis at x = 16 km (right) extracted from the 100 simulated velocity models.

Figure 8 .

 8 Figure 8. 68% and 95% isovalues of the marginal probability density function for the anisotropy parameters η and δ (left), P-velocity and S-velocity (middle), δ and the interface depth of reflector h5 (right). The represented values are the relative ones (for instance, η-η sol |ση| , where σ η is the standard deviation for the anisotropy parameter η). The crosses represent the MP values of the 100 simulated models.

Figure 9 .

 9 Figure 9. Solution model 2 obtained by tomography with additional constraints on the interface depth at well locations. The RMS value of the traveltime misfits is 6.0 ms. The anisotropy parameters values are η = 8.27% and δ = 15.94%. The crosses represent the wells.

Figure 10 .

 10 Figure 10. Solution model 3 obtained by tomography with a constrain on δ parameter to remain positive. The RMS value of the traveltime misfits is 6.371 ms. The anisotropy parameters values are η = 6.2% and δ = 2%.

Table 1 .

 1 Standard deviations in bold (square roots of the diagonal terms of the a posteriori covariance matrix) associated with the MP 'mean in the illuminated part of the model' and the normalized correlations between the uncertainties on the different defined MP. This uncertainty analysis is performed for the solution model of figure1. The normalized values

	σ |value sol |	× 100 of

Table 2 .

 2 RMS values of the traveltime misfits, mean value of the velocities VP5 and VS5, of the interfaces H 4 and H 5 in the illuminated region (MP defined in section 4.2 and anisotropy parameters η and δ for 20 simulated models (chosen among the 100 simulations presented in section 4.1).Models RMS (ms) VP5 (km s -1 ) VS5 (km s -1 ) H 4 (km) H 5 (km) η (%) δ (%)

	1	6.7	3.736	1.154	1.783	3.131	6.5 -2.6
	2	6.6	3.121	1.222	1.925	3.222	6.2 -4.9
	3	6.4	2.746	1.102	1.790	3.060	6.0 -5.9
	4	6.9	2.701	1.315	1.857	3.259	6.2 -3.7
	5	6.2	3.300	1.416	1.817	3.142	5.7 -8.3
	6	6.5	2.542	1.245	1.926	3.181	6.1 -4.6
	7	6.6	2.536	1.137	1.925	3.143	6.1 -5.5
	8	6.2	1.919	1.087	1.844	3.102	6.3 -4.8
	9	6.2	3.182	1.078	1.875	3.267	6.5 -3.5
	10	6.4	2.960	1.388	1.736	3.198	6.0 -6.2
	11	6.2	3.247	1.452	1.815	3.173	6.0 -6.4
	12	6.2	2.592	1.073	1.664	3.147	6.6 -1.4
	13	6.5	2.510	1.239	1.817	3.239	6.2 -4.9
	14	6.6	3.595	1.130	1.717	3.208	6.5 -1.9
	15	11.6	2.344	1.172	1.770	3.144	6.0 -6.1
	16	6.5	2.154	1.323	1.853	3.2322	6.3 -4.5
	17	6.3	2.855	1.307	1.797	3.296	6.3 -3.4
	18	6.2	2.125	1.049	1.921	3.118	5.9 -6.8
	19	6.6	1.431	1.048	1.772	3.144	6.4 -3.7
	20	6.3	2.371	1.225	1.864	3.166	6.0 -6.6

Table 3 .

 3 Computational times (platform: Silicon Graphics Fuel V10 Graphics 500 MHz R14000A) for unconstrained inversion, constrained inversion, Cholesky factorization and one macro-parameter (MP).

		Unconstrained	Constrained	Cholesky	Macro-parameter
		inversion	inversion	factorization	(MP)
	CPU	1797 s/iteration	1930 s/iteration	9251 s	1684 s/MP
	time	with 1779 s/iteration	with 1779 s/iteration	for 100 simulated	without
		for ray tracing	for ray tracing	models without	ray tracing
				raytracing

Table 4 .

 4 Three different solution models. Model 1, model 2 and model 3 are, respectively, the models of figures 1, 9 and 10.

		Model 1	Model 2	Model 3
	Mean of the P-velocity	2.849 km s -1	3.623 km s -1	2.895 km s -1
	Mean of the S-velocity	1.229 km s -1	1.621 km s -1	1.213 km s -1
	η	6.29%	8.27%	6.2%
	δ	-4.4%	15.94%	2%
	Mean of the depth of the reflector H 4	1.838 km	1.755 km	1.659 km
	Mean of the depth of the reflector H 5	3.167 km	3.021 km	2.882 km

In our approach (see Jurado et al[START_REF] Jurado | 3D reflection tomography designed for complex structures Expanded Abstracts[END_REF]), the subsurface model m is composed of 2D or 3D B-spline functions describing velocity variations in a layer (v(x, y) + k.z or v(x, y, z)) and 2D B-spline functions describing the interfaces (Z(x, y), Y (z, x) or X(y, z)) (see[START_REF] De Boor | A Practical Guide to Splines[END_REF][START_REF] Inoue | A least-squares smooth fitting for irregularly spaced data: finite element approach using the cubic B-splines basis[END_REF]).

The computation of the derivatives of T with respect to the model parameters is cheap, thanks to the Fermat principle (see[START_REF] Bishop | Tomographic determination of velocity and depth in laterally varying media[END_REF]). Thus, the derivatives of traveltimes with respect to velocity parameters are calculated by integrating a perturbation of the velocity along the rays calculated in the background model. The derivatives with respect to interface parameters are the integration of the background velocity along the shortened (negative contribution) or extended (positive contribution) part of the ray resulting from the interface displacement (we consider only the vertical displacement of impact point). See[16] for details.

PS-wave results from the conversion at the reflector of a down-going P-wave (compressional wave) into an up-going S-wave (shear wave). In opposition to the PS-wave, the pure P mode is often called the PP-wave.

The expressions of η and δ are δ= (C 13 +C 44 ) 2 -(C 33 -C 44 ) 2 2C 33 (C 33 -C 44 )and η = -δ 1+2δ where = C 11 -C 33 2C 33 , C ij being the elements of the elastic stiffness tensor C, which defines the linear relationships between the stress tensor and the strain tensor (following[START_REF] Thomsen | Weak elastic anisotropy[END_REF] and[START_REF] Stopin | Reflection tomography of compressional and shear modes for determination of anisotropic velocity models[END_REF]).

From figure 2, we define the illuminated region by x ∈ [14, 20.5 km].

Acknowledgments

The authors would like to thank BP for providing the data set, Karine Broto, Frédéric Delbos, Andréas Ehinger, Anne Jardin, Patrick Lailly and Timothée Perdrizet for many fruitful discussions.