Julien Gance
email: j.gance@brgm.fr

Gilles Grandjean

Kévin Samyn

Jean-Philippe Malet

Quasi-Newton inversion of seismic first arrivals using source finite bandwidth assumption: application to landslides characterization from velocity and attenuation fields

Keywords: Inverse theory, Seismic tomography, Fresnel wavepaths, Finite frequency, Resolution improvement, Landslide

Introduction

Near-surface geophysical techniques are of great potential to image subsurface structures in many geomorphological studies (Schrottand Sass, 2008;[START_REF] Van Dam | Landform characterization using geophysics-Recent advances, applications, and emerging tools[END_REF], and even more specifically for improving the knowledge on landslide features (JongmansandGarambois, 2007). Landslide analysis generally involves the combined use of several geophysical methods to image different petrophysical parameters [START_REF] Bruno | Test of high-resolution seismic reflection and other geophysical techniques on the Boup landslide in the Swiss Alps[END_REF][START_REF] Mauritsch | Geophysical investigations of large landslides in the Carnic Region of southner Austria[END_REF]Méric et al., 2011). Seismic surveys are particularly well adapted to detect and characterize specific features of landslides such as spatial variation in compaction and rheology of the material [START_REF] Grandjean | Geophysical data fusion by fuzzy logic for imaging the mechanical behavior of mudslides[END_REF], variation in fissure density at the surface [START_REF] Grandjean | Characterization of a landslide fissure pattern by integrating seismic azimuth tomography and geotechnical testing[END_REF], and complex slip surface geometries [START_REF] Grandjean | Caractérisation de la structure interne et de l'état hydrique de glissements argilo-marneux par tomographie géophysique :l'exemple du glissement-coulée de Super-Sauze (Alpes du Sud, France)[END_REF]. As the wave propagation is mainly controlled by the elastic properties of the medium, seismic surveys are often well correlated with geotechnical observations. Among them, different approaches can be used for processing data depending on the analysisof the different kind of waves associated to particular propagation phenomena. [START_REF] Bichler | Threedimensional mapping of a landslide using a multi-geophysical approach: the Quesnel Forks landslide[END_REF] and [START_REF] Bruno | Test of high-resolution seismic reflection and other geophysical techniques on the Boup landslide in the Swiss Alps[END_REF] interpreted the late arrivals of P-waves to image the bedrock geometry using seismic reflection analysis. Mauristch et al. (2000) and [START_REF] Glade | Determination of potential landslide shear plane depth using seismic refraction. A case study in Rheinhessen, Germany[END_REF] analyzed refracted waves to define the internal geometry (e.g. layering) of the landslides. [START_REF] Grandjean | Geophysical data fusion by fuzzy logic for imaging the mechanical behavior of mudslides[END_REF] and [START_REF] Grandjean | Structural study of the Ballandaz landslide (French Alps) using geophysical imagery[END_REF] studied the first arrival traveltimes to recover P-wave velocity distribution across landslides. More recently, Samyn et al. (2011) used a 3D seismic refraction traveltime tomography to provide a valuable and continuous representation of the 3D structure and geology of a landslide. [START_REF] Grandjean | Geophysical data fusion by fuzzy logic for imaging the mechanical behavior of mudslides[END_REF] and [START_REF] Grandjean | Structural study of the Ballandaz landslide (French Alps) using geophysical imagery[END_REF] used also Spectral Analysis of Surface Waves (SASW) techniques to obtain S-wave velocity distributions along landslide cross-sections.

Finally, based on the work of [START_REF] Pratt | Seismic waveform inversion in the frequency domain. Part 1: Theory and verification in a physic scale model[END_REF] and VirieuxandOperto (2009), [START_REF] Romdhane | Shallowstructure characterization by 2D elastic full-waveform inversion[END_REF] demonstrated the possibility to exploit the entire wave signal by performing a Full elastic Waveform Inversion (FWI) on a dataset acquired on a clayey landslide.

These methods are based on more or less important assumptions on the soil type and their petrophysical characteristics, on the complexity of layers, and on their geometry, and are therefore not recommended for all case studies. Seismic refraction can be basically performed using the General Reciprocal Method (GRM; [START_REF] Palmer | The generalized reciprocal method of seismic refraction interpretation[END_REF] or the Plus/Minus Method [START_REF] Hagedoorn | The plus-minus method of interpreting seismic refraction sections[END_REF]. These methods are fast and easy to implement, but mostly based on simple hypotheses -such as a layered media -and are consequently not very efficient to recover important variations in lateral velocity. Reflection surveys need to record high frequency waves to obtain a good resolution of final stacked sections; they are difficult to set up for subsurface applications because of the strong attenuation affecting high frequencies and decreasing the signal to noise ratio (Jongmansand Garambois, 2007). Moreover, this method produces a reflection image that cannot be directly and easily interpreted asdirect information on the mechanical properties of the medium that could be derived from velocity fields. FWI appears to be the more advanced method since it is based on the realistic assumption of an elastic media and uses the whole seismic signal in the inversion scheme. However, it remains a complex method, requiring an important data pre-processing (source inversion, amplitude correction, initial model estimation) that is hardly applicable to near surface real data featured with a low signal to noise ratio and complex waves interactions [START_REF] Romdhane | Shallowstructure characterization by 2D elastic full-waveform inversion[END_REF]. As a consequence, the issue of recovering the structural image of a landslide from the seismic velocity field estimated with an accurate method is still a challenge.

The objective of this work is to refine the first arrival tomography approach applied to landslide analysis, which is a good compromise between the strong assumptions made in simple refraction methods and the complexity of the FWI approach when used in very heterogeneous soils. P-wave tomography is largely used, from crustal to shallow geotechnical studies and allows recovering, at a relatively low cost, a reliable image of the velocity distribution. In this perspective, Taillandier et al. (2009) proposed to invert first arrival traveltime tomography using the adjoint state method, but had to face the problem of gradient regularization. The proposed method is based on a Hessian formulation [START_REF] Tarantola | Inverse problem theory.Methods for data fitting and model parameters estimation[END_REF] to ensure an optimum convergence of the velocity model during iterations. We only use here the first arrivals of the seismic signal caused by direct or refracted waves. In that case, late arrival, including scattering, conversion effect are not taken into account in the inversion.

In the next section, we detail the theory used to formulate the inversion equation and the related approximations introduced for optimizing the algorithm and ensuring the convergence toward the best solution. The final algorithm is tested on a synthetic case in order to estimate the quality of the velocity field reconstruction without considering noisy data. To compare our results with a reference case, the synthetic velocity model is identical to the one used by [START_REF] Romdhane | Shallowstructure characterization by 2D elastic full-waveform inversion[END_REF]. Then, we process a real seismic dataset obtained from a seismic survey at the Super-Sauze landslide (South French Alps) in clayey material. On this case study, we invert for the P-wave velocity field but also for the attenuation factor. The effects of noisy data on the quality of the reconstructed velocity field are analyzed, and the observed surface variation in the velocity field is discussed by integrating other sources of information (geotechnical tests and geomorphological interpretation).

Theoretical approach a. General inverse problem

Solving the inverse problem requires a modeling step (e.g. the direct problem) for computing the residuals (e.g. the difference between computed and observed data). An updated model can be estimated by back-projecting the data residuals on each cell of the model discrete grid. The direct problem thus consists in finding a relation between the data t, taken here as the first arrival traveltimes, and the physical P-wave velocity models, such as:

𝑡 = 𝑓(𝑠) (1)
This problem is generally solved using the asymptotic high frequency wave propagation assumption in a non-homogeneous isotropic medium. The wave propagation equation is then simplified into the eikonal equation that computes the first arrival traveltimes over a discrete grid for example by using a finite-difference scheme (Vidale, 1998). This technique has been improved for sharp velocity contrasts (PodvinandLecomte, 1991), optimized by PopoviciandSethian (1998) who proposed the Fast Marching Method (FMM), by [START_REF] Zhao | A fast sweeping method for eikonal equations[END_REF] who proposed the Fast Sweeping Method, and recently adapted to non-uniform grids [START_REF] Sun | A finite difference scheme for solving the eikonal equation including surface topography[END_REF]. The technique is actually largely used because of its efficiency and accuracy. In this study, we use the FMM algorithm (PopoviciandSethian, 1998) to compute the direct problem. The computation of the first arrival traveltime is generally followed by a ray tracing that uses the traveltimes maps to propagate the wavepath from the sources to the receivers. The traveltime t can then be expressed using the following matrix form:

𝑡 = 𝐿𝑠 (2)
where L is the length of the ray segment crossing each cell of the slowness model S. Generally, the Fermat principle allows the Fréchet derivatives matrix to be approximated by the matrix 𝐿 [START_REF] Baina | Tomographie sismique entre puits : mise en oeuvre et rôle de l'analyse à posteriori vers une prise en compte de la bande passante[END_REF]. This principle establishes that a small perturbation 𝛿𝑠 on the initial model does not change the ray paths, and therefore, only affect the traveltime 𝛿𝑡 at the second order. This linearization of the problem around the slowness model 𝑠 𝑘 implements the Quasi-Newton (Q-N) assumption and makes the resolution of the nonlinear problem to be solved by the following tomographic linear system:

𝐿 𝑘 𝑇 𝐿 𝑘 𝛿𝑠 = 𝐿 𝑘 𝑇 𝛿𝑡 (3)
where k and T represent respectively the iteration index and the transpose operator. Several techniques can be used to solve Eq. 3. The Simultaneous Iterative Reconstruction Technique (SIRT) is the most common because it does not require a large computer memory; consequently, SIRT is widely used to invert large sparse linear systems (Trampertand Leveque, 1990). Its convergence toward a least-square solution has been proved (Van der Sluisand [START_REF] Van Der Sluis | Numerical solutions of large, sparse linear systems arising from tomographic problems[END_REF] but this method suffers from few drawbacks: SIRT introduces an intrinsic renormalization of the tomographic linear system that modifies slightly the final solution. For that reason, we propose to use in our approach the LSQR (Least-Square QR) algorithm [START_REF] Paige | LSQR: An algorithm for sparse linear equations and sparse least squares[END_REF] that has been proved to be superior to SIRT or Algebraic Reconstruction Techniques ART algorithms in terms of numerical stability and convergence rates [START_REF] Nolet | Solving or resolving inadequate and noisy tomographic systems[END_REF].

b. From rays to Fresnel volumes

Revisiting seismic traveltimes tomography also needs to consider the ray approximation, e.g. the infinite spectral bandwidth assumption. The main issue of this approximation lies in the traveltime computation taken as line integrals along the rays spreading over the slowness model. Because the slowness values are only considered along ray paths, the problem is often underdetermined and leads to numerical instability [START_REF] Baina | Tomographie sismique entre puits : mise en oeuvre et rôle de l'analyse à posteriori vers une prise en compte de la bande passante[END_REF]. In practice, this difficulty is generally by-passed using regularization operators to reduce the non-constrained part of the model, and then, to reinforce the numerical stability. Generally, a simple smoothing of the reconstructed slowness model (Zeltand Barton, 1998), the application of a low-pass filter on the gradient (Taillandier, 2009) or the elimination of the lowest eigenvalues in the Hessian matrix [START_REF] Tarantola | Inverse problem theory.Methods for data fitting and model parameters estimation[END_REF]) can be used. However, such regularization operators require the selection of appropriate parameters (filter length, eigenvalue cutoff and more generally the size and weights of the smoothing operators). [START_REF] Kissling | Model parameterization in seismic tomography: a choice of consequence for the solution quality[END_REF] demonstrated the dependency of these parameters on the resolution and the quality of the final model in poor wavepath coverage areas. This is the reason why the use of physically-based regularization operators such as Fresnel volumes or sensitivity kernels is preferable since they are based on nonsubjective principles completely defined by the problem.

With this assumption, several methods have been developed including those using the concept of Fresnel volumes as a regularization factor. [START_REF] Nolet | Seismic wave propagation and seismic tomography[END_REF] proposed to use the size of the Fresnel zone to constrain the size of a spatial smoothing operator. [START_REF] Vasco | Beyond ray tomography: Wavepaths and Fresnel volumes[END_REF], [START_REF] Watanabe | Seismic traveltime tomography using Fresnel volume approach[END_REF], [START_REF] Grandjean | JaTS: a fully portable seismic tomography software based on Fresnel wavepaths and a probabilistic reconstruction approach[END_REF] or [START_REF] Ceverny | Seismic ray theory[END_REF] use the Fresnel volume in the back projection of the residuals.

More recently, the concept of sensitivity kernels or Fréchet Kernel introduced by [START_REF] Tarantola | Inverse problem theory.Methods for data fitting and model parameters estimation[END_REF] has been reformulated for traveltime tomography. It represents a good compromise between the strong assumption of the asymptotic ray theory and the high computational cost of the complete full wave inversion tomography [START_REF] Liu | Sensitivity kernels for seismic Fresnel volume tomography[END_REF]. While ray theory is well adapted in media characterized by geological structures larger than the first Fresnel zone, the use of sensitivity kernels allows to overcome this constraint and then, to increase the spatial resolution (SpetzlerandSnierder, 2004).

Several authors have developed this concept and investigated the properties of the sensitivity kernels for traveltime tomography in homogeneous media, mainly at the spatial scale of the earth crust [START_REF] Dahlen | Fréchet kernels for finite-frequency traveltimes-I[END_REF][START_REF] Dahlen | Finite-frequency sensitivity kernels for boundary topography perturbations[END_REF][START_REF] Zhao | Structural sensitivities of finite-frequency seismic waves: A fullwave approach[END_REF]2006). The strategies proposed by these authors are generally limited to smoothly heterogeneous media [START_REF] Liu | Sensitivity kernels for seismic Fresnel volume tomography[END_REF][START_REF] Spetzler | Case story: Timelapse seismic crosswell monitoring of CO2 injected in an onshore sandstone aquifer[END_REF] and are barely applied to real datasets, mainly because of prohibitive calculation costs [START_REF] Liu | Sensitivity kernels for seismic Fresnel volume tomography[END_REF].

Sensitivity kernels and Fresnel volumes are constructed on the single scattering (Born) approximation. This approximation considers that a part of the wave can be delayed by velocity perturbation in the vicinity of the ray path. For highly heterogeneous media, this approximation is not valid anymore because multiple scattering should be taken into account to fully explain observed first arrivals.

In the particular case of subsurface soil investigation where highly heterogeneous media are observed, the use of sensitivity kernels is not recommended since it would involve complex algorithms and computing times similar to those featuring the FWI approach. In this context, we decided to solve empirically the problem of complex multiple scattering inside the first Fresnel zone, this approach allowing to address issues of regularization by using increasing finite frequencies bandwidths.

Therefore, the developed algorithm assumes thathigh spatial heterogeneities of the soil around the Fresnel wavepath affect, through complex multiple scattering, the first arrivals of the signal.

In the next section, we explain how the Fresnel volumes are defined as a simplification of the sensitivity kernels and implemented in a Q-N algorithm.

We first express the Fresnel weights proposed by [START_REF] Watanabe | Seismic traveltime tomography using Fresnel volume approach[END_REF] and used by [START_REF] Grandjean | JaTS: a fully portable seismic tomography software based on Fresnel wavepaths and a probabilistic reconstruction approach[END_REF], in order to materialize the wavepath in the model. They classically decrease linearly from a value of one (when the cell is positioned on the ray path) to a value of zero (when it is out of the Fresnel volume): [START_REF] Watanabe | Seismic traveltime tomography using Fresnel volume approach[END_REF] and [START_REF] Grandjean | JaTS: a fully portable seismic tomography software based on Fresnel wavepaths and a probabilistic reconstruction approach[END_REF] defined these weightsω as the probability that a slowness perturbation delay the arrival of the wave by a Δt, where f is the considered frequency of the wave. This definition is interesting because of its simple expression that allows fast computing while taking into account the global shape of the 2D traveltime sensitivity kernel given of SpetzlerandSnieder ( 2004), corresponding to a decrease of sensitivity until the first Fresnel zone.

𝜔 = 1 -2𝑓Δ𝑡, 0 ≤ Δ𝑡 < 1 2𝑓 0, 1 2𝑓 ≤ Δ𝑡 (4)
Because the size of the Fresnel volume is depending on the considered source frequency, a new inversion strategy based on increasing frequency can be proposed. To consider the entire source frequency content, some authors compute wavepaths in a band-limited sensitivity kernel, stacking the monochromatic kernels with a weight function similar to the amplitude spectrum of the wavelet (SpetzlerandSnieder, 2004;[START_REF] Liu | Sensitivity kernels for seismic Fresnel volume tomography[END_REF]. From our side, we propose to compute the Fresnel weights for a monochromatic wave, increasing its frequency at each step of the inversion. The considered frequency sampling with an increasing rule ranging from the lower to the higher frequency of the source signal. We also choose to give the same weight to all frequencies in order to limit the number of iteration and to preserve the rapid convergence of the algorithm. Taillandier et al. (2009) showed that this method, applied on the gradient filtering permit to overcome the non-linearity of the problem. For low frequency values, the Fresnel zone will be wider and the slowness model will be reconstructed with large wavelengths. Conversely, for high frequency values, the Fresnel zone will be thinnerand the slowness model will be reconstructed with sharp wavelengths. This strategy appears efficient to improve the resolution during the inversion -and so the convergence of the algorithm -in full wave inversion where it prevents from cycle-skipping issues (Sirgueand Pratt, 2004;VirieuxandOperto, 2009;[START_REF] Romdhane | Shallowstructure characterization by 2D elastic full-waveform inversion[END_REF]. Recovering slowness variations whose sizes are in agreement with each wavelength transmitted by the source should lead to better images of local heterogeneities while preserving the algorithm convergence.

With these assumptions, the traveltime perturbation can be expressed as the integral over a Fresnel volume multiplied by the slowness perturbation field observed for all points r in the volume [START_REF] Liu | Sensitivity kernels for seismic Fresnel volume tomography[END_REF]. This linear relationship among traveltime and slowness perturbation is the result of the first Born approximation that is only valid for small perturbation [START_REF] Yomogida | Fresnel zone inversion for lateral heterogeneities in the Earth[END_REF]:

𝛿𝑡 = 𝑊 𝑟 𝛿𝑠 𝑟 𝑑𝑟 (5)
where𝑊 is the Fresnel weight operator (Fig. 1) normalized by the area of the Fresnel ellipse perpendicular to the raypath𝑎, for each shot and receiver such as:

𝑊 𝑟 = 1 𝑎 𝜔(𝑟) (6)
This normalization allows linking our approach to the geometrical ray theory in the case of a plane wave propagation [START_REF] Vasco | Beyond ray tomography: Wavepaths and Fresnel volumes[END_REF]Spetzler and Snieder, 2004;[START_REF] Liu | Sensitivity kernels for seismic Fresnel volume tomography[END_REF], so that the integral of W over a surface perpendicular to the wavepath is equal to one, which is generally verified for sensitivity kernels:

𝑊 𝑠, 𝑟, 𝑓 𝑑𝑠 = 1 +∞ -∞ (7)
The Fresnel volume thus defined, also called Fréchet kernel corresponds to the Fréchet derivative [START_REF] Yomogida | Fresnel zone inversion for lateral heterogeneities in the Earth[END_REF][START_REF] Tarantola | Inverse problem theory.Methods for data fitting and model parameters estimation[END_REF] or Jacobian matrix of the forward problem, similar to the length of ray segment matrix 𝐿 defined in the asymptotic ray theory. The problem, considered in 2D, can then be written in its matrix form [START_REF] Vasco | Beyond ray tomography: Wavepaths and Fresnel volumes[END_REF]:

𝛿𝑡 = 𝑊 𝛿𝑠 (8)
where𝑊 is the [number of datax number of cells] Fresnel weight matrix, calculated for each shot and for each receiver such as for the cell j and for the couple source-receiver i:

𝑊 𝑖𝑗 = 𝜔 𝑖𝑗 𝑙 (9)
where𝑙 is the length of the Fresnel surface along a direction perpendicular to the ray path.

We are however aware that the use of the Fresnel weight as defined previously is a rough approximation of the wavepath. This is the reason why we prefer to keep the ray assumption to compute traveltimes, avoiding the introduction of any error due to this approximation. Thischoice has been justified by [START_REF] Liu | Sensitivity kernels for seismic Fresnel volume tomography[END_REF] who observed small differences between traveltime values computed with both approaches.

c. New implementation of the inverse Problem

In the previous section, we showed that the Born approximation combined with Fresnel volumes allowed the estimation of Fréchet derivative. To solve the nonlinear inverse problem, we use a steepest descent iterative algorithm to minimize the 𝐿 2 norm misfit function l [START_REF] Tarantola | Inverse problem theory.Methods for data fitting and model parameters estimation[END_REF] that can be written in this particular case:

𝑙 𝑠 = 1 2 𝑓(𝑠) -𝑡 𝑜𝑏𝑠 𝑇 𝐶 𝑇 -1 𝑓(𝑠) -𝑡 𝑜𝑏𝑠 + 𝑠 -𝑠 𝑝𝑟𝑖𝑜𝑟 𝑇 𝐶 𝑆 -1 (𝑠 -𝑠 𝑝𝑟𝑖𝑜𝑟 ) (10)
where𝐶 𝑇 and 𝐶 𝑆 are respectively the covariance operators on data and model, 𝑓 represents the theoretical relationship between the model 𝑠 and the traveltimes 𝑡, 𝑡 𝑜𝑏𝑠 the observed data and 𝑠 𝑝𝑟𝑖𝑜𝑟 the a priori information on the model.

The Q-N method consists in minimizing the misfit function iteratively using its gradient and approximated Hessian matrix. The gradient gives the direction of steepest descent of the misfit function while the Hessian matrix is used as a metric indicating its curvature, approximated locally by a paraboloïd [START_REF] Tarantola | Inverse problem theory.Methods for data fitting and model parameters estimation[END_REF]. Compared to the gradient method, the Hessian is here used to improve the assessment of direction and norm of the update vector that is applied to the slowness model, so that each step is performed in the good direction and with an optimized length. The tomographic linear system in the case of no a priori information on the model can be written in its matrix form:

𝑊 𝑘 𝑇 𝐶 𝑇 -1 𝑊 𝑘 𝛿𝑠 = 𝑊 𝑘 𝑇 𝛿𝑡 (11)
This linear tomographic system can be solved by using a LSQR algorithm. The Fresnel weights matrix employed is a [number of data x number of cells] matrix. Its computation is straightforward on classical personal computers. The Q-N algorithm is based on the inverse of the Hessian matrix that is a square matrix of size [number of data x number of data] that can be difficult to invertfor large dataset.

In our case, we used cluster computing to invert our data and optimized the step size along the direction given by the Hessian using a scalar to weight the slowness update thanks to a parabolic interpolation [START_REF] Tarantola | Inverse problem theory.Methods for data fitting and model parameters estimation[END_REF].

d. Validation on a synthetic dataset

To compare the performance of our algorithm with full wave inversion, it was tested on the synthetic transverse section used by [START_REF] Romdhane | Shallowstructure characterization by 2D elastic full-waveform inversion[END_REF] representing a typical cross-section of a landslide (Fig. 2a). The performance of the method is discussed by taking as a reference the SIRT algorithm of [START_REF] Grandjean | JaTS: a fully portable seismic tomography software based on Fresnel wavepaths and a probabilistic reconstruction approach[END_REF].

The synthetic dataset has been calculated with a simple 2D eikonal equation solver, so that traveltimes are not perturbed by scattering effects. The dataset is composed of 50 shots recorded on 100 geophones. The model consists in 209x68 cells of 1 m in width. The seismic sources are spaced regularly every 2 m and the geophones every 1 m. The synthetic model is composed of ten layers with P-wave velocity ranging from 110 m.s -1 to 3300 m.s -1 . It represents a transversal cross-section of the landslide body lying on a homogeneous consolidated bedrock. Except between the landslide body and the bedrock, the changes in P-wave velocity are small, so that it is difficult to reconstruct the detailed shape of the landslide body. The stopping criterion is a change in the cost function lower than 1% and the inversion is limited to 20 steps.

Figure 2b and 2c compare the results obtained with the SIRT algorithm [START_REF] Grandjean | JaTS: a fully portable seismic tomography software based on Fresnel wavepaths and a probabilistic reconstruction approach[END_REF] and the Q-N one. They converge respectively in 15 and 20 iterations. The convergence of the misfit function is better for the Q-N algorithm with a value 3.6 times lower than the SIRT. Globally, the Q-N algorithm succeeds in recovering the shape of the bump at 2000 m.s -1 located on the left side of the cross-section and renders a more spatially detailed shape of the bedrock topography (Fig. 2b,2c). The SIRT algorithm does not retrieve those initial structures. We can notice that both algorithms have difficulties to reproduce the low velocity of the very near surface layer, because of the poor ray coverage occurring in this area.

Figure 2d, 2e and 2f compare the results on three vertical cross-sections chosen to highlight the ability of the algorithm to recover P-wave velocity with a better accuracy than the SIRT. We can notice that for the top-soil, P-wave velocity are the same for both algorithms and are equal to the P-wave velocities of the initial model in that zone. Moreover, The Q-N algorithm seems to be able to recover lowest wavelengths than the SIRT one as on figure 2d between 55 and 60 m of depth.

This synthetic example demonstrates the ability of the Q-N algorithm to recover the shape of the internal layers of a landslide body. Contrary to the SIRT result, the Q-N allows to interpret correctly the thickness and geometry of the different layers.

Application of the Q-N algorithm on a real dataset a. Study site presentation

The Q-N algorithm has been applied on a real dataset acquired at the Super-Sauze landslide (South French Alps). The landslide has developed in Callovo-Oxfordian black marls. Its elevation is between 2105 m at the crown which is established in in-situ black marls covered by moraine deposits and 1740 m at the toe. The landslide is continuously active with displacement rates of 0.05 to 0.20 m.day -1 (Malet et al., 2005a). The detachment of large blocks of marls from the main scarp (Fig. 3) and their progressive mechanical and chemical weathering in fine particles explain the strong grain size variability of the material especially in the topsoil (e.g. decametric blocks of marls at various stages of weathering, decimetric and centimetricclasts of marls, silty-clayed matrix; [START_REF] Maquaire | Instability conditions of marlyhillslopes: towards landsliding of gullying? The case of the Barcelonnette Basin, South East France[END_REF]. The bedrock geometry is complex with the presence, in depth, of a series of in-situ black marl crests and gullies, partially or totally filled with the landslide material [START_REF] Flageollet | The 3D Structure of the Super-Sauze Earthflow: A first stage towards Modelling its Behavior[END_REF][START_REF] Travelletti | Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data[END_REF]. This complex bedrock geometry delimits sliding compartments of different hydrogeological, rheological and kinematical pathways. The variable displacement rates and the bedrock geometry also control the presence of large fissures at the surface (Fig. 3b,3c) that can be imaged with joint electrical and electromagnetic methods [START_REF] Schmutz | Joint electrical and time domain electromagnetism (TDEM) data inversion applied to the Super-Sauze earthflow (France)[END_REF]. Many kinds of heterogeneities are observed at different scales, and they control directly the mechanical behavior of the landslide by creating excess pore water pressures [START_REF] Van Asch | Influence of landslide geometry and kinematic information to describe the liquefaction of landslides: Some theoretical considerations[END_REF][START_REF] Travelletti | Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data[END_REF]. The top soil surface characteristics have also a large influence on the surface hydraulic conductivity [START_REF] Malet | Soil surface characteristics influence on infiltration in black marls: application to the Super-Sauze earthflow (Southern Alps, France)[END_REF] and therefore, on the water infiltration processes (Malet et al., 2005b;[START_REF] Debieche | Local scale groundwater modelling in a landslide. The case of the Super-Sauze mudslide (Alpes-de-Haute-Provence, France)[END_REF]. In this context, we tested the Q-N P-wave tomography inversion scheme to provide a high-resolution characterization of structures and detect small scale heterogeneities along aN-S transect of the landslide.

b. P-wave velocities inversion

The seismic profile is parallel to the main sliding direction of the material, and is located in the upper part from the secondary scarp to the middle part of the accumulation zone (Fig. 3). The base seismic device has a length of 94 m long and consisted of 48vertical geophones (with a central frequency of 10 Hz) regularly spaced every 2 m. The 60 shots have been achieved with a hammer every 4 m. We used a roll-along system to translate the acquisition cables with a 24 geophones overlap and then investigate a 238 m long linear profile. The recording length is 1.5 s with a sample rate of 0.25 ms and the acquisition central consists in a Geometrics Stratavizor seismic camera (48 channels). The shots show a relatively good signal to noise ratio until the last geophone in the upper part of the section. The shots of the lower part are affected by the higher density of soil fissures. The attenuation and the signal to noise ration are consequently higher in this part of the profile. The signal is dominated by the surface waves but the first arrivals are clearly visible for all seismic shots.

To estimate the pick uncertainty, the timedifference of the picks from reciprocal source-receiver pairs was examined. Two different people picked the first arrival of the waves, and for each created dataset we studied the difference between reciprocal traveltimes. From each dataset, we created 3 subsets: a first one kept intact, a second one where we corrected the picks having more than 10 ms of difference in traveltime reciprocity, the picks being removed when not possible and a third one removing picks having a difference of traveltime reciprocity greater than 10 ms. We inverted the 6 datasets and kept the one which gave the lower final misfit function value and the P-wave tomogram the most reliable with our a priori knowledge of the studied zone. The differences between P-wave tomograms were not so important, but the dataset selected permitted to interpret the results unequivocally. The error in reciprocity of traveltimes can be approximated by a Gaussian belt with a standard deviation of 2.2 ms (Fig. 5).

The inversion has been performed with the Q-N algorithm on a 358 x 153 grid. Each square cell of the grid measures 0.67 m. The hammer source gave frequencies comprised between 30 Hz and 120 Hz with a dominant frequency of 40 Hz. Those frequencies are also present in the P-wave up to the end of the profile (Fig. 4b). Regarding the different amplitude spectra, we considered those values as constant for each shot. For the first iteration, the frequency was set to 30 Hz and then increases at each step to respectively 45 Hz, 60 Hz, 75 Hz, 90 Hz, 105 Hz and 120 Hz; after the eighth steps, the frequency was kept constant at 120 Hz.

The initial model is of a first importance for a good convergence. Indeed, the linearization of the problem around the initial model is only valid for model close to the real one. This reason motivates the automated calculation of the initial model from the data by using a simplified slant-stack algorithm transforming the dataset in the velocity -intercept times domain. Then, theoretical traveltime curves were calculated in a double loop for a velocity v rangingfrom 𝑣 𝑚𝑖𝑛 to 𝑣 𝑚𝑎𝑥 (defined by the user) and for a range of intercepts 𝜏 such as:

𝑡 = 𝜏 + 𝑜𝑓𝑓𝑠𝑒𝑡 𝑣 (12)
For each of these theoretical curves, the model was compared to the observed traveltimes; the number of points with a difference of less than 0.5 ms was stored in the 𝜏-p plan (Fig. 6a). For each intercept, the velocity corresponding to the maximum number of point was used to create a profile of velocity that is further used under the shot point (Fig. 6b).

The depth 𝑑 corresponding to each velocity was then deduced from the velocities by using Eq. 13:

𝑑 𝑖 = 𝑑 𝑖 -1 + 𝑣 𝑖 * 𝛿𝑡 (13)
where𝛿𝑡 represents the sample rate of intercepts (Fig. 6c).

The velocity profiles were gathered to form the initial model that was finally smoothed to avoid lateral artifacts. With a misfit function value of 0.3612, this initial model is supposed to be as much as possible close to the final one which is the best way to avoid divergence or convergence in local minimum. The stopping criterion was set to a percentage change lower than 1 % and the number of iteration was limited to 20. The misfit function value associated to the final model is 0.0122 and has been obtained after 17 iterations.Figure 7 shows the misfit functions related to the Q-N and SIRT algorithms for comparison.

Figure 8 shows respectively the initial velocity model and the inverse models obtained with both SIRT and Q-N algorithms.

c. Seismic wave attenuation tomography

The amplitude of the first arrival is exploited here to image the wave attenuation. Seismic wave attenuation is an efficient physical property because it is directly linked to porosity and to the presence of fissures in the media [START_REF] Schön | The physical properties of rocks: fundamentals and principles of petrophysics[END_REF]. Several attenuation parameters can be used to invert seismic attenuation: spectral ratio, centroid frequency shift and peak frequency shift [START_REF] Castro Nunes | Estimating quality factor from surface seismic data: a comparison of current approaches[END_REF]. In our analysis, we used the simple method proposed by [START_REF] Watanabe | Seismic attenuation Tomography and its application to rock mass evaluation[END_REF] considering that in a homogeneous attenuation media the amplitude of the spherical wave verifies:

𝐴 𝑟 = 𝐴 0 𝑟 𝑒 -𝛼𝑟 (14)
where𝐴 0 is the amplitude of the source signal, r the distance from the seismic source along the raypath and α the attenuation coefficient.

Applying a logarithm function, the equation becomes linear in 𝛼. We proposed to compute the attenuation tomography after the P-wave velocity field so that the problem is reduced to a simple linear problem:

𝛼 = 𝛼 + 𝑊 𝑇 𝛿𝐴 ( 15
)
where α is the attenuation coefficient, 𝛿𝐴, the amplitude update and W the Fresnel weights matrix.

We performed five iterations starting from a simple homogeneous media with an attenuation 𝛼 = 1.0 𝑒 -03 Np.m -1 and calculated the Fresnel weights matrix for the Q-N inverted Vp model and for five frequencies (30 Hz,45 Hz,60 Hz,90 Hz,120 Hz). Although the results are generally presented through the dimensionless quality factor, in our case, the dominant source frequency is almost constant for each shot (40 Hz) and the raw results are more contrasted. Moreover, because the definition of the quality factor can be different for lossy (highly attenuating) materials [START_REF] Schön | The physical properties of rocks: fundamentals and principles of petrophysics[END_REF], we image directly the attenuation map. In order to compare the attenuation map with surface fissuring, a geomorphological inventory of the fissures was created by direct observations in the field along the seismic profile. All the fissures of width larger than 0.05 m were mapped. The Surface Cracking Index (SCI), defined as the total length of fissures per linear meter of the seismic profile, is represented in Figure 9a.

Interpretation and discussion

The velocity model is firstly compared to the geometrical model proposed by [START_REF] Travelletti | Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data[END_REF] proposed from the integration of multi-source data (e.g. electrical resistivity tomographies, dynamic penetration tests and geomorphological observations) at coarser scale.

The bedrock depths obtained from penetration tests were used to deduce the minimum velocity of the bedrock fixed at 850 m.s -1 . The bedrock position thus obtained is compared to the one of [START_REF] Travelletti | Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data[END_REF] along the profile (Fig. 10a,b). The thickness was calculated for both models with the topography observed in 2011. The depth interval observed is the same (2 to 13 m) and the global shape of the bedrock is identical. The main differences are lower than 4 m and are located downhill of the profile, from the abscissa 100 to 230 m where the bedrock topography is more complex. It shows important high frequencies variations that are less pronounced on the topography from [START_REF] Travelletti | Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data[END_REF].

To assess the quality of the Q-N algorithm result, we compare the bedrock depth variation with the topography visible on a black and white orthophotograph of 1956 before the landslide (e.g. when the original relief was not covered by the landslide material). The figure 11a shows the orthophotograph of 1956 overlaid on the digital elevation model available for the same date. The seismic profile is also represented. We can notice on the picture the very irregular geometry of the relief composed of alternating crests and gullies of different sizes. The studied profile crosses different geometries. From the point A to B (Fig. 11a), the topography seems first regular, before crossing a large thalweg. After crossing this thalweg, the profile goes across a long transversal crest that forms a bump and then, plunges along the downstream slope. Those patterns are also present in the P-wave tomogram, respectively at 125, 140 m and 160 m.

Near the abscissa 180 m, a zone of low velocities (< 900 m.s -1 ) is observed in depth that can be interpreted as an old small-sized gully, also visible as a black dot on the orthophotograph.

From this observation, we constructed an interpreted geological model that gathers geophysical and geomorphological information (Fig. 11b). This model is composed of three different materials. The first one is the Callovo-Oxfordian intact black marls that constitute the bedrock, characterized by velocities greater than 3000 m.s -1 [START_REF] Grandjean | Geophysical data fusion by fuzzy logic for imaging the mechanical behavior of mudslides[END_REF] and with a topography following the one visible on the DEM from 1956. The second one is the unconsolidated material that presents P-wave velocities between 400 and 1200 m.s -1 , characteristic of unconsolidated deposits [START_REF] Bell | Engineering Geology[END_REF]. It constitutes the active unit of the landslide, e.g. the layer that presents most displacements. The third material, characterized by an interval of velocities between 1500 and 2500 m.s -1 is interpreted as compacted landslide material, because it is only present at a higher depth, where the height of soil above is sufficient to compact the landslide deposit. This layer is not involved in the dynamics of the landslide andis quasi-impermeable (Malet, 2003;[START_REF] Flageollet | The 3D Structure of the Super-Sauze Earthflow: A first stage towards Modelling its Behavior[END_REF][START_REF] Travelletti | Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data[END_REF].

Now that we proved the reliability of the Q-N algorithm tomography to recover the paleotopography, we can assess its contribution by comparing this result from the SIRT tomography of [START_REF] Grandjean | JaTS: a fully portable seismic tomography software based on Fresnel wavepaths and a probabilistic reconstruction approach[END_REF]. We can notice on figure 8b and 8c, that the upper part of the profile (between the abscissa 0 and 100 m) presents a smoother contrast between the active unit and the bedrock for the SIRT result. Although two bumps are visible at the abscissa 100 and 140 m, the large in-between depression is not visible. In the lower part, the active layer is deeper than for the Q-N algorithm tomography and presents little lateral heterogeneity that prevent the interpreter from delimitating the landslide layer with accuracy. Globally, the SIRT tomography and the Q-N one are in agreement, but SIRT algorithm does not permit to recover the geotechnical unit limits with the same accuracy and prevents from distinguishing the "Dead body" represented in figure 11b.

To summarize, we proved that the Q-N algorithm previously developed presented a better horizontal and vertical resolution than the SIRT algorithm [START_REF] Grandjean | JaTS: a fully portable seismic tomography software based on Fresnel wavepaths and a probabilistic reconstruction approach[END_REF] and gives the possibility to interpret the bedrock geometry as geomorphological and geological patterns with more details and better resolution than previously. Thanks to the algorithm, we were able to clearly distinguish three of the four geotechnical units proposed by [START_REF] Flageollet | The 3D Structure of the Super-Sauze Earthflow: A first stage towards Modelling its Behavior[END_REF]. The active unit composed of C1a and C1b geotechnical layers, the "dead body" that represents the C2 layer and the bedrock.

The seismic attenuation field obtained is in agreement with the P-wave velocities inversed. Although P-wave velocity and seismic attenuation are not linked during the inversion, the attenuation field also permits to distinguish to main materials. The first one presents a relatively low attenuation, that correspond to the bedrock and the "dead body" gathered. The values obtained for this material are around 10 -3 Np.m -1 , which corresponds to consolidated soils attenuation at 40 Hz [START_REF] Schön | The physical properties of rocks: fundamentals and principles of petrophysics[END_REF]. The second one is a more attenuating and heterogeneous material that corresponds to the landslide deposit.

The attenuation of the landslide layer varies from 4 to 12 10 -3 Np.m -1 which are more typical values for unconsolidated soils at this frequency [START_REF] Schön | The physical properties of rocks: fundamentals and principles of petrophysics[END_REF]. The lateral heterogeneity of the attenuation is important. Areas of high attenuation show attenuation values three times higher than others. The correlation between the SCI and the attenuation is very good, so that we can affirm that attenuation variations are mainly caused by cracking.

The seismic wave attenuation tomography highlights two important attenuation zones. The first one, in the upper-part at the abscissas 20-30 seems to be concentrated at the interface between the landslide layer and the bedrock (Fig. 11b). In this area the landslide material undergoes significant shear stresses due to friction between the landslide layer and the bedrock topography. At the surface, cracks wider than 15 cm and measuring several meters in length are visible and the soil is completely remolded.The second zone seems to be associated to the influence of the local bedrock geometry that forms a bump followed by a steep slope (Fig. 11b). Its concave shape is responsible of the tension cracks well visible at the surface and clearly perpendicular to the major change of topography direction (Fig. 3c).

Moreover, the soil seems dryer in that zone, certainly because of the geometry of the area that facilitates the drainage of the groundwater table toward lower areas. This phenomenon can even intensify the presence of cracks.

We can notice that P-wave velocity and seismic wave attenuation produce complementary results.

Only the P-wave velocity tomography can reproduce the sharp geometry of the bedrock topography and differentiate the "dead body" from the bedrock (Fig. 11b). On the other hand, the seismic wave attenuation provides more information on the landslide layer lateral heterogeneity.

Here, we can notice that the two high attenuation zones are not correlatedwith low P-wave velocity as expected [START_REF] Grandjean | Structural study of the Ballandaz landslide (French Alps) using geophysical imagery[END_REF]. This is explained by the seismic wavepaths going through this weathered material, that are probably too short to impact the velocity, so that only scattering and diffusion effects affect the wave amplitudes.

Those complementary results highlight the importance of seismic attenuation tomography when characterizing weathering state of soils in near surface application.

It is finally important to note that this geophysical method is, in theory, not adapted to image the near surface zone (first meters in depth). Indeed, the poor Fresnel weighs coverage in that zone prevents to obtain high accuracy results. Moreover, cracks and attenuating objects size is certainly comparable or smaller than the cell size (0.67 m). So, we have to be cautious when interpreting those results because the real attenuation zones could be closer to the surface than those presented on the seismic attenuation tomography.

Conclusions

We present a P-wave tomography inversion algorithm for highly heterogeneous media which combine improvement of resolution and inversion regularization. We use the Fresnel wavepaths calculated for different sources frequencies to retropropagate the traveltime residuals, assuming that in highly heterogeneous media, the first arrivals are affected by velocity anomalies present in the first Fresnel zone through complex multiple scattering. After verifying the ability of the algorithm to recover complex structures on a synthetic dataset, we apply it on a real dataset acquired on the upper part of the Super-Sauze landslide. We verify that our results are in accordance with previous results Such geological models are generally created for hazard assessment through numerical modeling.

For this work, the use of geophysical tomography is essential because it is the only tool that can provide continuous and integrated imaging of the soil. In this perspective, our algorithm contribution is to answer one of the major issues raised by [START_REF] Travelletti | Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data[END_REF], namely, the lack of resolution of geophysical tomography compared to geotechnical tests and geomorphological observations. This is why we believe that the Q-N algorithm could really improve the resolution of the geological model, and thus any numerical hydro-mechanical modeling using this variable as a parameter. Interpreted geological cross-section in 3 layers: the active unit corresponds to the moving landslide layer, the dead body is a compacted and quasi-impermeable layer showing low displacements and the 754 bedrock constituted of in-situ black marls.

Figure Captions

from

  [START_REF] Travelletti | Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data[END_REF] and dynamics penetrometer tests. It appears that the sharp geometry of the bedrock topography is well recovered. This result is of first importance because bedrock topography is one of the main controlling factors of landslide displacement. Using the wavepaths calculated for P-wave velocity tomography inversion, we invert separately the seismic wave attenuation of the same dataset. Although barely applied in landslide environment, we show that with few efforts it was possible to use the amplitude of the first peak of the wave to get an accurate image of the seismic wave attenuation of the landslide layer. The results are in agreement with the observed surface crack inventory and bring complementary information for the construction of an interpreted model. The use of traveltimes and amplitude of waves permit to create an interpreted geological model that gathers those different details.

Figure 1 :

 1 Figure 1: Fresnel weights (e.g. distance versus depth) computed for a medium characterized by a

Figure 2 :

 2 Figure 2: Algorithm validation on a synthetic dataset: a) Synthetic initial model, b) Final velocity

Figure 3 :

 3 Figure 3: Location of the investigated area within the Super-Sauze landslide. The seismic profile

Figure 4 :

 4 Figure 4: Data quality: a) Seismic shot# 13. b) Source amplitude spectrum of the total seismic shot

Figure 5 :

 5 Figure 5: Picking quality: Frequency distribution of differential traveltime errors (e.g. due to

Figure 6 :

 6 Figure 6: Example of vertical P-wave velocity profile for the seismic shot #4: a) Slant-stack

Figure 7 :Figure 8 :

 78 Figure 7: Misfit function values for the Q-N and SIRT algorithms. 731

Figure 9 :

 9 Figure 9: Seismic wave attenuation tomography: a) Surface Craking Index (SCI) and b)

Figure 10 :

 10 Figure 10: Comparison of the bedrock geometry interpreted with the Q-N algorithm and with a

Figure 11

 11 Figure 11: a) Pre-event topography before the landslide and location of the seismic profile, b)

Acknowledgements

This work was supported by the French National Research Agency (ANR) within the Project SISCA "'SystèmeIntégré de Surveillance de Crises de glissements de terrains Argileux, 2009-2012" and the BRGM Carnot institute.


G Blekherman

J B Lasserre
email: lasserre@laas.fr

The truncated K-moment problem for closure of open sets *

Keywords: Moment problem, truncated moment problem, maximum entropy. Subject class: 44A60 65K10 42B10

We consider the truncated K-moment problem when K ⊆ R n is the closure of a, not necessarily bounded, open set (which includes the important cases K = R n and K = R n + ). We completely characterize the interior of the convex cone of finite sequences that have a representing measure on K ⊆ R n . It is in fact the domain of the Legendre-Fenchel transform associated with a certain convex function. And so in this context, detecting whether a sequence is in the interior of this cone reduces to solving a finite-dimensional convex optimization problem. This latter problem is related to maximum entropy methods for approximating an unknown density from knowing only finitely many of its moments. Interestingly, the proposed approach is essentially geometric and of independent interest, as it also addresses the abstract problem of characterizing the interior of a convex cone C which is the conical hull of a set continuously parametrized by a compact set M ⊂ R n or M ⊆ S n-1 , where M is the closure of an open subset of R n (resp. S n-1 ). As a by-product we also obtain a barrier function for the cone C.

Introduction

We are concerned with the (real) truncated K-moment problem, that is, given a closed set K ⊆ R n , and a finite sequence y = (y α ), α ∈ N n 2d (where N n 2d = {α ∈ N n : i α i ≤ 2d}), provide conditions under which y has a representing Borel measure on K, i.e., y is such that

y α = K x α dµ(x), ∀α ∈ N n 2d ,
for some finite Borel measure µ on K.

Background: For the one-dimensional (or univariate) case n = 1, this classical problem is well understood and dates back to contributions by famous mathematicians, among them Markov, Stieltjes, Hausdorff, and Hamburger, at the end of the nineteen and beginning of twentieth centuries. Explicit conditions on the sequence y exist, all stated in terms of positive semidefiniteness of some Hankel matrices whose entries are linear in the variables y (see e.g. Curto and Fialkow [START_REF] Curto | Recursiveness, positivity and truncated moment problems[END_REF]); in modern language, these conditions are Linear Matrix Inequalities (in short LMIs) in y.

For the multi-dimensional (n > 1) case, no such strong results exist, even for the full moment problem, and for instance, the full moment problem with K = R n is still unsolved. The Riesz-Haviland criterion states that an infinite sequence y = (y α ), α ∈ N n , has a representing measure on K if and only if y (viewed as a linear functional acting on the polynomials) is nonnegative for all polynomials nonnegative on K; but since there is no tractable characterization of the latter polynomials, the Riesz-Haviland criterion is not practical. Existence of a representing measure is related to existence of commuting self-adjoint extensions of (multiplication) operators on polynomials, defined from the sequence y (see e.g. Berg [START_REF] Berg | The multidimensional moment problem and semi-groups[END_REF], Sarason [START_REF] Sarason | Moment problems and operators in Hilbert space[END_REF], Simon [START_REF] Simon | The classical moment problem as a self-adjoint finite difference operator[END_REF], Vasilescu [START_REF] Vasilescu | Spectral measures and moment problems, Spectral Theory and Its Applications[END_REF]) and so far, the most powerful (and general) result is due to Schmüdgen [START_REF] Schmüdgen | The K-moment problem for compact semi-algebraic sets[END_REF], who solved the full K-moment problem when K is a compact basic semi-algebraic set of the form K := {x ∈ R n : g j (x) ≥ 0, j = 1, . . . , m} for some polynomials (g j ) ⊂ R [x]. In this context, a sequence y has a representing measure on K if and only if it satisfies countably many (explicit) LMI's; this result was later refined (and simplified) by Putinar [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF] when the quadratic module generated by the g j 's is Archimedean. Later, the full K-moment problem for basic closed (not necessarily compact) semialgebraic sets was also solved (at the price of a dimensional extension) in Putinar and Vasilescu [START_REF] Putinar | Solving the moment problem by dimensional extension[END_REF]. Finally, there also exist conditions in terms of linear inequalities on y, based on an alternative representation theorem initially due to Krivine [START_REF] Krivine | Anneaux préordonnés[END_REF][START_REF] Krivine | Quelques propriétés des préordres dans les anneaux commutatifs unitaires[END_REF], and also later in Marshall [START_REF] Marshall | A general representation theorem for partially ordered commutative rings[END_REF] and Vasilescu [START_REF] Vasilescu | Spectral measures and moment problems, Spectral Theory and Its Applications[END_REF].

However, for the truncated moment problem in a general context, the only "explicit" criterion is the so-called flat extension of positive moment matrices in Curto and Fialkow [START_REF] Curto | The truncated complex K-moment problem[END_REF][START_REF] Curto | Truncated K-moment problems in several variables[END_REF]. (A positive semidefinite moment matrix M d (y) associated with y ∈ N n 2d has a flat extension if the sequence y can be extended to ỹ ∈ N n 2d+2 in such a manner that the resulting moment matrix M d+1 (ỹ) has same rank as M d (y).) Namely, y = (y α ), α ∈ N n 2d , has a representing measure on R n if M d (y) is positive semidefinite and y can be extended to ỹ ∈ N n 2(d+k) for some k, in such a manner that ỹ has a flat extension. But again this test is not practical. Finally, in Jordan and Wainwright [START_REF] Wainwright | Graphical Models, Exponential Families, and Variational Inference[END_REF] the authors characterize the sequences y that have a representing measure with a density with respect to a reference measure.

Contribution: We consider the truncated K-moment problem where K ⊆ R n is the closure of a not necessarily bounded open set; and so, in particular, it includes the important special cases K = R n and K = R n + . In this context, we completely characterize the interior of the convex cone

C(K) ⊂ R s(n,2d) (where s(n, d) := n+d n ) of finite sequences y = (y α ), α ∈ N n
2d , that have a finite representing measure on K. Namely, let µ be any measure on K, absolutely continuous with respect to the Lebesgue measure on K. Then, if K is compact we show that any sequence y ∈ int(C(K)) has a representing measure ν absolutely continuous with respect to µ, and such that

y α = K x α e p(s) dµ(s) dν(s) , ∀α ∈ N n 2d ,
for some polynomial p ∈ R[x] of degree at most 2d. We note that the above criterion depends only on the moments of degree at most 2d, and does not rely on extending the moment matrix or taking higher moments into account.

If K = R n then we take a little detour in R n+1 by homogenization, so that µ is now a rotation invariant measure on the unit sphere S n ⊂ R n+1 , and y ∈ int(C(K)) ⊂ R s(n+1,2d) has a representing measure ν absolutely continuous with respect to µ, and such that

y α = S n
x α e p(s) dµ(s)

dν(s) , ∀α ∈ N n+1 2d , (1.1)
for some homogeneous polynomial p ∈ R[x 0 , • • • , x n ] of degree 2d. Importantly, the above result holds, with identical proof via homogenization, when K is the closure of an unbounded open subset of R n (like e.g. R n + ). In this case, and as explained and detailed in Section 3.3, integration on S n in (1.1) is replaced with integration on some compact subset M of S n .

Alternatively, y ∈ int(C(K)) has a representing measure if and only if f * (y) < +∞, where f * : R s(n,2d) → R ∪ {+∞} is the Legendre-Fenchel transform of the convex function

p → f (p) := K e p(s) dµ(s) or M e p(s) dµ(s) , (1.2)
defined for polynomials of degree at most 2d. That is, f * is defined as

y → f * (y) := sup p∈R[x] 2d {p T y -f (p) }.
And so, checking whether y ∈ int (C(K)) reduces to solving the finitedimensional convex optimization problem P of finding the supremum of p T y -f (p). We show that the supremum is finite and attained on the interior of the cone C(K) and f * (y) = +∞ for y not in the interior. This means that f * (resp. log f * ) provides a barrier (resp. log-barrier) function for the cone C(K).

Our result is in the vein of (and extends) Wainwright and Jordan [START_REF] Wainwright | Graphical Models, Exponential Families, and Variational Inference[END_REF]Theorem 3.3], where (when the domain of f is open) the authors have shown that the gradient map ∇ log f is onto the interior of the convex set M of sequences y that have total mass 1 and a representing measure with a density with respect to a reference measure µ; here we prove the same result for the interior of the (potentially larger) convex cone of sequences that have a representing measure not necessarily absolutely continuous w.r.t. µ. Furthermore, we guarantee that the Legendre-Fenchel transform f * is finite only on the interior of the cone of representable moment sequences, while [START_REF] Wainwright | Graphical Models, Exponential Families, and Variational Inference[END_REF]Theorem 3.3] makes no guarantees of the behavior of (log f ) * on the boundary. This extra regularity of f * allows us to conclude that it is a barrier function for the cone C(K).

However, and even though P is a finite dimensional convex problem, effective numerical computation of f * (y) is still difficult. This is because evaluating f and its gradient ∇f (as well as its Hessian ∇ 2 f for secondorder methods) at a point y, requires evaluating integrals over K, a difficult problem. However, notice that if K is relatively "simple", one may approximate those integrals by using cubature formulas, or discretization schemes, or Monte-Carlo methods; also, for small dimensions 2 and 3, and if K is defined by polynomials, then these integrals can be approximated efficiently and accurately by dimension reduction to line and surface integrals; see e.g. Wester et al. [START_REF] Wester | Computing integrals over polynomially defined regions and their boundaries in 2 and 3 dimensions[END_REF].

The optimization problem P is well-known and called the maximumentropy approach (in our case, the Boltzmann-Shannon entropy) for estimating an unknown density from the only knowledge of finitely many of its moments. In maximum entropy the sequence y is known to come from a representing measure and the main question of interest is the convergence of an optimal solution p * d ∈ R[x] 2d of P when the number of moments (i.e., d) increases. For a detailed account of such results, the interested reader is referred to Borwein and Lewis [START_REF] Borwein | On the convergence of moment problems[END_REF] and the many references therein. And so, another contribution of this paper is to show that the maximum entropy approach not only permits to approximate an unknown density but also permits to solve the K-moment problem for closure of open sets.

Finally, our (essentially geometric) approach is also of independent interest, as we obtain our result on the truncated K-moment problem as a by-product of the more abstract problem of characterizing the interior of a convex cone C ⊂ R m which is the conical hull of a compact set L(M) continuously parametrized by a set M ⊂ R n or M ⊆ S n , where M is the closure of an open subset of R n (resp. S n-1 ).

Notation, definitions and preliminary results

Let V be the Euclidean space with an inner product , and L : R n → V a continuous mapping. Given a compact set M ⊂ R n that is the closure of a bounded open subset of R n or S n-1 , let L(M) ⊂ V be a compact subset lying in an affine hyperplane H ⊂ V .

Let h ∈ V be the vector perpendicular to H such that x, h = 1 for all x ∈ H. Let µ be a measure on M absolutely continuous with respect to the Lebesque measure or the rotation-invariant measure on S n-1 , and with a density positive on M, and normalize µ to have mass 1.

We are primarily interested in the conical hull of L(M), which is the convex cone

C := ConicalHull(L(M)) = λ i x i | λ i ≥ 0, x i ∈ L(M) .
Without loss of generality we may assume that the affine hull of L(M) is all of H, so that the cone C is full-dimensional in V . Define the function f : V → R as follows:

x → f (x) := M e x,L(v) dµ(v), x ∈ V. (2.1)
The function f is smooth and strictly convex (follows from the definition of µ and M compact). Using the compactness of M we may differentiate under the integral sign to obtain:

∇f (x) = M L(v) e x,L(v) dµ(v), x ∈ V.
It follows that f is a function of "Legendre-type" in the sense of Rockafellar [START_REF] Rockafellar | Convex Analysis[END_REF]Chapter 26,p. 258]. If we think of ∇f as a function mapping V to V , then it is clear that the image of ∇f lies in C. Define the Legendre-Fenchel transform f * : V → R of f as follows:

y → f * (y) := sup x∈V x, y -f (x), y ∈ V.
By [13, Theorem 26.5], ∇f is one-to-one and the image of ∇f is the interior of the domain of f * (i.e., points where f * is finite). In particular, it follows that the image of ∇f is convex. Moreover, the inverse of ∇f , viewed as a mapping from V to V , is just ∇f * , i.e.:

(∇f ) -1 = ∇f * .

Next, we introduce a couple of intermediate lemmas that we will need later to prove our main result. They follow from elementary convexity and analysis and for clarity of exposition, their proofs are postponed until Section 4.

Lemma 2.1. Let B be a compact convex set and let A be a convex subset of B such that the closure A of A contains all exposed extreme points of B. Then A contains the interior of B. Lemma 2.2. Let M be a compact set that is the closure of an open subset of R n or S n-1 , and let µ be a measure on M, absolutely continuous with respect to the Lebesque or rotation invariant measure on S n-1 , and with a density positive on M. Suppose that f : M → R is a continuous function such that f is positive on M and it attains its maximum at a unique point s ∈ M. Then for all continuous functions g : M → R,

lim λ→∞ M gf λ dµ M f λ dµ = g(s).
(2.2)

We now show that the image of ∇f is the interior of C.

Theorem 2.3. The image of ∇f is the interior of C.

Proof. From the above it follows that the image of ∇f is an open convex set. Observe that for all a ∈ R

f (x + ah) = M e x+ah,L(v) dµ(v) = e a M e x,L(v) dµ(v), ∀x ∈ V. Therefore, ∇f (x + ah) = e a ∇f (x), ∀a ∈ R, ∀x ∈ V.
It follows that the image of ∇f is an open convex sub-cone of C.

Next, let B = conv L(M) (the convex hull of L(M)) be the base of the cone C. It suffices to show that ∇f is onto the interior of B. By Lemma 2.1 it is enough to show that we can approximate any exposed extreme point of B arbitrarily well by points ∇f (x i ) for some sequence {x i } ⊂ V . So, let s ∈ B be an exposed extreme point. It follows that s ∈ L(M) and we can write s = L(s ′ ) with s ′ ∈ M. Furthermore there exists p ∈ V such that p, L(s ′ ) = 0 and p, L(v

) < 0 for all v ∈ M with v = s ′ . Now consider the point ∇f (β p), β ∈ R + , i.e., ∇f (β p) = M L(v) e β p,L(v) dµ(v), β ∈ R + .
To make sure that the points we consider lie in B we need to divide by M e β p,L(v) dµ(v). Define

s β := M L(v) e β p,L(v) dµ(v) M e β p,L(v) dµ(v)
.

Applying Lemma 2.2 with f := e p,L(v) yields lim β→∞ s β = s. Therefore, the sequence of points ∇f (β p+a h), β ∈ R + , with e -a := M e β p,L(v) dµ(v), approximate s, the desired result.

We know from [START_REF] Rockafellar | Convex Analysis[END_REF] that the image of ∇f is the interior of the domain of f * . We have shown that on the interior of the image of ∇f the supremum is always attained. Also, outside of C we know that f * is equal to +∞. Now we show that under an additional assumption on the geometry of the embedding L(M), f * = +∞ on the boundary of C. We will show that the additional assumption holds in our applications of interest. α -e α f (β p).

Consider f (β p) = M e β p,L(v) dµ. We know that e p,L(v) is at most 1 on M and by the assumption of the Lemma, the maximum of 1 is attained on a set of measure zero. Therefore f (βp) can be made arbitrarily small by taking β appropriately large.

It follows that m α,β = +∞, since we can take arbitrarily large α and then adjust β so that e α f (βp) is arbitrarily small. Thus f * (y) = +∞ for all y ∈ ∂B.

Moment Cones

Now we apply the geometric machinery we developed to the moment cones.

The (compact) truncated K-moment problem.

Let V := R[x] 2d denote the vector space of polynomials in n variables of degree at most 2d, and let K ⊂ R n be a compact set that is the closure of an open subset of R n . For this application the set M of Section 2 will simply be K. Let µ be a measure, supported on K, that is absolutely continuous with respect to the restriction of the Lebesque measure to K.

Let P 2d (K) denote the cone of polynomials in R[x] 2d that are nonnegative on K. The cone

P 2d (K) is a closed, convex, full-dimensional, pointed cone in R[x] 2d . Let R 2d (K) ⊂ R[x] *
2d denote the cone of linear functionals on R[x] 2d that come from integration with respect to a finite Borel measure supported on K, i.e. the set of all linear functionals ℓ ∈ R[x] * 2d which can be written in the form

ℓ(p) = K p dσ
for some measure σ.

For every v ∈ R n , let ℓ v ∈ R[x] * 2d denote the linear functional given by evaluation at the point v: ℓ v (p) = p(v). We can view ℓ v as the integrational functional with respect to the Dirac-δ measure on v. Let L(K) ⊂ R 2d (K) denote the set of all linear functionals

ℓ v ∈ R[x] * 2d with v ∈ K. Clearly, the set L(K) is a continuous embedding of K into R[x] * 2d taking v ∈ K to ℓ v . The function f in (2.1) now reads p → f (p) := K e ℓv,p dµ(v) = K e p(v) dµ(v), p ∈ R[x] 2d .
Let H ⊂ R[x] * 2d be the affine hyperplane of all linear functionals that evaluate to 1 on the constant polynomial 1:

H = {ℓ ∈ R[x] * 2d | ℓ(1) = 1}.
It is clear that L(K) is contained in H. Finally, in order to apply our framework we claim that R 2d (K) is the conical hull of L(K).

Lemma 3.1. The cone R 2d (K) of all linear functionals representable by a measure supported on K is the conical hull of the set L(K) of linear functionals ℓ v with v ∈ K.

Proof. Let C(K) denote the conical hull of L(K). We note that since L(K) is compact and included in the hyperplane H it follows that the cone C(K) is closed. The inclusion C(K) ⊆ R 2d (K) is straightforward. To prove the reverse inclusion note that the dual cone C(K) * of C(K) is the cone P 2d (K) of polynomials non-negative on K, simply because the functionals ℓ v with v ∈ K encode non-negativity on K. By bi-duality it follows that C(K) = P * 2d (K). For every ℓ ∈ R 2d (K), ℓ(p) ≥ 0 for all p ∈ P 2d (K) and therefore R 2d (K) ⊆ P * 2d (K). Thus we obtain

C(K) ⊆ R 2d (K) ⊆ P * 2d (K) = C(K),
which yields the desired result.

Lemma 3.1 can also be derived from extensions of Tchakaloff's theorem in Putinar [START_REF] Putinar | A note on Tchakaloff's theorem[END_REF] and Bayer and Teichman [START_REF] Bayer | The proof of Tchakaloff's theorem[END_REF]. Now we can directly apply Theorem 2.3 to R 2d (K) with M = K. It still remains to check that the assumptions of Lemma 2.4 hold for L(K). Let F be a maximal (by inclusion) face of R 2d (K). Maximal faces of a convex cone are exposed. Therefore, there exists a form p ∈ P 2d (K) such that ℓ(p) = 0 for all ℓ ∈ F and ℓ(p) > 0 for all ℓ ∈ R 2d (K) \ F . Now suppose that L(v) = ℓ v ∈ F . It follows that ℓ v (p) = p(v) = 0. Therefore the set of v ∈ K for which L(v) is in F corresponds precisely to the zeroes of p in K:

{v ∈ K | L(v) ∈ F } = {v ∈ K | p(v) = 0}.
Since K is the closure of an open set in R n and the measure µ is absolutely continuous with respect to the Lebesque measure it follows that µ({v

∈ K | L(v) ∈ F }) = 0.
In summary we have proved the following result:

Theorem 3.2. Let K ⊂ R n be the closure of an open bounded subset. Let µ be an arbitrary finite Borel measure on K, absolutely continuous with respect to the restriction of the Lebesgue measure on K, with a density positive on K, and let

p → f (p) := K e p(x) dµ(x), p ∈ R[x] 2d . (3.1) A sequence y ∈ R[x] * 2d belongs to int R 2d (K) if and only if f * (y) := sup p { p, y -f (p) } < +∞. (3.2)
In other words, int R 2d (K) is the domain of the Legendre-Fenchel transform of f .

Observe that the function f * (resp. log f * ) provides a barrier (resp. log-barrier) for the convex cone R 2d (K).

Interestingly, the function log f is well-known to statisticians. It is called the log partition (or, cumulant) function associated with the so-called potential functions (or, sufficient statistics) (x α ), α ∈ N n 2d . Rephrased in our context, Theorem 3.3 in Wainwright and Jordan [START_REF] Wainwright | Graphical Models, Exponential Families, and Variational Inference[END_REF] states that under weak hypotheses, the mapping ∇(log f ) is onto the open convex set of moment sequences y ∈ R[x] * 2d that have a representing measure ν absolutely continuous with respect to µ and total mass 1. Our result is an extension of [START_REF] Wainwright | Graphical Models, Exponential Families, and Variational Inference[END_REF]Theorem 3.3] as we prove that ∇f is onto int R 2d (K), i.e., the interior of the cone of sequences that have arbitrary representing measures (as opposed to measures absolutely continuous w.r.t. µ). Furthermore, we guarantee that the Legendre-Fenchel transform f * is finite only on the interior of the cone of representable moment sequences, while [START_REF] Wainwright | Graphical Models, Exponential Families, and Variational Inference[END_REF]Theorem 3.3] makes no guarantees of the behavior of (log f ) * on the boundary.

The truncated moment problem on R n .

We now consider the more delicate case when K = R n . In this case, as K is not compact the above machinery cannot be applied directly with M = K and a detour is needed.

Let R 2d (R n ) ⊂ R[x] * 2d denote the cone of linear functionals on R[x] 2d that come from integration with respect to a finite, Borel measure supported on R n , i.e. the set of all linear functionals ℓ ∈ R[x] * 2d which can be written in the form

ℓ(p) = R n
p dσ for some measure σ.

Unlike the situation of the compact support K, the cone R 2d (R n ) is no longer closed. However, there is a nice way to represent the closure of R 2d (R n ). It is well-known that the dual cone of R 2d (R n ) is the cone P 2d (R n ) of polynomials nonnegative on all of R n . By bi-duality, the closure of R 2d (R n ) is the cone P * 2d (R n ), i.e.:

P * 2d (R n ) = R 2d (R n ).
Given an arbitrary polynomial p ∈ R[x] 2d we can homogenize p by adding an extra variable x 0 and multiplying all monomials in p by an appropriate power of x 0 so that all monomials have degree 2d. Let p denote the homogenization of p. Conversely, we can de-homogenize p by setting x 0 = 1 to obtain p. If p is a nonnegative polynomial, then p is a nonnegative form.

Let V = H n,2d denote the vector space of all homogeneous forms in n + 1 variables of degree 2d. We can linearly identify R[x] 2d with H n,2d via homogenization. Let HP 2d denote the cone of nonnegative forms on R n+1 . From the above it follows that homogenization identifies the cone P 2d (R n ) with the cone HP 2d .

Define HR 2d ⊂ H * n,2d to be the cone of linear functionals on H n,2d given by integration with respect to a finite Borel measure on R n+1 . For this example the compact set M of Section 2 will be the unit sphere S n in R n+1 . For v ∈ R n+1 let ℓ v ∈ H * n,2d be the linear functional given by evaluation at v, i.e.:

p → ℓ v (p) = p(v) for all p ∈ H n,2d .
As before, let L(M) be the set of linear functionals ℓ v with v ∈ M (= S n ). Let H be the hyperplane in H * n,2d consisting of all functionals that evaluate to 1 on (x 2 0 + . . . We next show in the following Lemma that HR 2d is a closed convex cone and in fact HR 2d is the conical hull of L(M). This is very similar to the situation in Lemma 3.1 and the proof is almost identical.

+ x 2 n ) d . It is clear that L(M)
Lemma 3.3. The cone HR 2d of all linear functionals representable by a measure supported on R n is the conical hull of the set L(M) of linear functionals ℓ v with v ∈ M = S n .

Proof. Let C(M) denote the conical hull of L(M). We note that since L(M) is compact and included in the hyperplane H it follows that the cone C(M) is closed. The inclusion C(M) ⊆ HR 2d is straightforward. To prove the reverse inclusion note that the dual cone C(M) * of C(M) is the cone HP 2d of non-negative forms, simply because the functionals ℓ v with v ∈ S n encode non-negativity on S n and by homogeneity on all of R n+1 . By bi-duality it follows that C(M) = P * 2d (M). For every ℓ ∈ HR 2d , ℓ(p) ≥ 0 for all p ∈ HP 2d and therefore HR 2d ⊆ HP * 2d . Thus we obtain C(M) ⊆ HR 2d ⊆ HP * 2d = C(M), which yields the desired result.

As we have seen in the proof of Lemma 3.3, the cone HR 2d is dual to the cone HP 2d of nonnegative forms. Via de-homogenization we can identify HP 2d with the cone of nonnegative polynomials P 2d (R n ). Therefore it follows that

HR 2d = R 2d (R n ) = P * 2d (R n ) = HP *
2d . Now we can apply Theorem 2.3 to the cone HR 2d . As before, we need to check that the assumptions of Lemma 2.4 hold for L(M). Let F be a maximal (by inclusion) face of HR 2d . Maximal faces of a convex cone are exposed. Therefore, there exists a form p ∈ HP 2d such that ℓ(p) = 0 for all ℓ ∈ F and ℓ(p) > 0 for all ℓ ∈ HR 2d \ F . Now suppose that L(v) = ℓ v ∈ F . It follows that ℓ v (p) = p(v) = 0. Therefore the set of v ∈ S n for which L(v) is in F corresponds precisely to the zeroes of p in S n :

{v ∈ S n | L(v) ∈ F } = {v ∈ S n | p(v) = 0}.
Since the measure µ is the rotation invariant probability measure on S n it follows that µ({v ∈ S n | L(v) ∈ F }) = 0.

In summary we have proved the following Theorem:

= y α , ∀(α, k) ∈ N n+1 2d , |α| + k = 2d,
has a representing measure on R n+1 . Indeed, y has a representing measure µ on R n if and only if ỹ has the representing (product) measure µ ⊗ δ x 0 =1 on R n+1 (where δ x 0 =1 is the Dirac measure at x 0 = 1).

General Non-Compact Case

With identical proofs via homogenization, the above discussion of the case K = R n can be extended to any set K which is the closure an open subset of R n . We explain how to define the appropriate compact set M ⊂ S n . Embed K into R n+1 by introducing an extra coordinate and setting it equal to 1. More formally, let K ′ ⊂ R n+1 given by:

K ′ = (x, 1) ∈ R n+1 | x ∈ K .
Next, define, M ′ ⊂ S n as a rescaling of K ′ onto the unit sphere:

M ′ = z ∈ S n | z = λy for some λ ∈ R, y ∈ K ′ .
Finally, let M ⊂ S n be the closure of M ′ : M = M ′ . Take any finite measure µ on M that is absolutely continuous with respect to the rotation invariant probability measure on the unit sphere S n . This makes sense, since K is the closure of an open subset of R n , and therefore M is a closure of an open subset of S n . We note that with K = R n the above construction gives M = S n . At last, nonnegativity of polynomials on K is equivalent to nonnegativity of forms on M. The rest of the proofs follow nearly word for word.

Proofs of Lemmas

Proof of Lemma 2.1. We observe that A ia a closed compact set. Since exposed extreme points are dense in the set of extreme points, it follows that A contains all extreme points of B and by Krein-Milman Theorem the condition of the lemma is equivalent to A = B. Now suppose that there exists x ∈ int B such that x / ∈ A. Then by the Separation Theorem, there exists a hyperplane H such that x lies in the closed half-space H + and A lies in the closed half-space H -.

Since x ∈ int B it follows that there exist points y ∈ B, such that y lies in the open half-space H -. Since A ⊂ H + we see that A cannot contain such points y, which is a contradiction.

Proof of Lemma 2.2. First we note that by dividing through by the maximum of f we may restrict ourselves to the case 0 < f (x) ≤ 1 for all x ∈ M. Next, for every λ ∈ N, define the probability measure g dν λ k , may be used to prove -ǫ ′ g(s) ≤ ρ -g(s) ≤ ǫ ′ g(s). Letting ǫ → 0 (hence ǫ ′ → 0) yields the desired result (2.2).

Lemma 2 . 4 .

 24 Suppose that for any face F of C we have µ({v | L(v) ∈ F }) = 0. Then f * (y) = +∞ for all y in the boundary of C. Proof. It suffices that show that f * (y) = +∞ for all y in the boundary ∂B of the base B. Let y ∈ ∂B. Then there exists p ∈ V such that p, y = 0 and p, L(v) ≤ 0 for all v ∈ M. Consider m α,β = sup α,β αh + βp, y -f (αh + βp). It follows that m α,β = sup α,β

  is a continuous embedding of S n into H * n,2d and L(M) lies in H. The analogue of the function f in (2.1) now reads p → f (p) = S n e ℓv,p dµ(v) = S n e p(v) dµ(v), p ∈ H n,2d .

Fix ǫ > 0

 0 and let M(ǫ) ⊂ M be the open set {x ∈ M : f (x) < 1 -ǫ}. As µ has a density positive on M, lim λ→∞ M(ǫ)f λ dµ 1/λ = esssup {f (x) : x ∈ M(ǫ)} = (1 -ǫ) (4M(ǫ)) 1/λ ,which in turn implies lim λ→∞ ν λ (M(ǫ)) = 0 for any fixed ǫ > 0. Therefore we also have lim λ→∞ ν λ (M \ M(ǫ)) = 1. Let us evaluateρ := lim sup λ→∞ M g dν λ ; ρ := lim inf λ→∞ M g dν λ .Let (λ j ), j ∈ N, be a subsequence such thatρ := lim j→∞ M g dν λ j , )| ν λ j (M(ǫ)), ∀ǫ > 0, so that lim j→∞ A ǫ j = 0. Therefore, j ) ν λ j (M \ M(ǫ)),for some ξ j ∈ M \ M(ǫ). Since lim λ→∞ ν λ (M \ M(ǫ)) = 1 we see thatρ = lim j→∞ g(ξ j ).Next, we also have(1 -ǫ ′ )g(s) ≤ g(x) ≤ (1 + ǫ ′ )g(s), ∀x ∈ M \ M(ǫ),with ǫ ′ → 0 as ǫ → 0. And so |g(ξ j ) -g(s)| ≤ ǫ ′ g(s) which implies-ǫ ′ g(s) ≤ ρ -g(s) ≤ ǫ ′ g(s).Exactly same arguments for the subsequence (λ k ), k ∈ N, such thatρ = lim inf λ→∞ M g dν λ = lim k→∞ M

  Theorem 3.4. Let µ be the rotation invariant probability measure on the unit sphere S n ⊂ R n+1 , and let In other words, int HR 2d is the domain of the Legendre-Fenchel transform of f .Observe that the function f * (resp. log f * ) provides a barrier (resp. logbarrier) for the convex cone HR 2d .Finally, to relate the initial moment problem in R n with the one in R n+1 with homogenization, observe that a sequence y = (y α ) ∈ R[x 1 , . . . , x n ] * 2d has a representing measure on R n if and only if the sequence ỹ ∈ H *

	p → f (p) :=	S n	e p(x) dµ(x),	p ∈ H n,2d .	(3.3)
	A sequence y ∈ H * n,2d belongs to int HR 2d if and only if
	f * (y) := sup p	{ p, y -f (p) } < +∞.	(3.4)
					n,2d
	defined by:
	ỹαk

Acknowledgement

The first author was partially supported by NSF and wishes to acknowledge financial support from the University of Toulouse, France, for a visit to LAAS in June-July 2011, during which this research was performed.