Modeste Irakarama
email: <modeste.irakarama@univ-lorraine.fr>

Paul Cupillard
email: 1<paul.cupillard@univ-lorraine.fr>

Guillaume Caumon
email: 1<guillaume.caumon@univ-lorraine.fr>

Paul Sava

Jonathan Edwards
email: <jonathan.edwards@univ-lorraine.fr>

Appraising structural interpretations using seismic

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

Uncertainties in structural interpretations of seismic images have been of interest to geoscientists for a long time. Earlier studies were mainly focused on migration velocity errors and seismic resolution; some examples can be found in [START_REF] Hajnal | Maximum uncertainty of interval velocity estimates[END_REF] and [START_REF] Beylkin | Spatial resolution of migration algorithms[END_REF]. Structural uncertainties propagated from velocity errors usually translate to subtle displacement or flexing of horizons and faults away from the reference position of those horizons and faults [START_REF] Thore | A practical formulation of migration errors due velocity uncertainties[END_REF][START_REF] Bube | The influence of stacking velocity uncertainties on structural uncertainties[END_REF][START_REF] Pon | Sensitivity analysis of seismic depth migrations[END_REF][START_REF] Fomel | Structural uncertainty of time-migrated seismic images[END_REF]. By the reference position of a given horizon or fault, one often means the position obtained from the tomographic velocity model (e.g. [START_REF] Messud | Extracting seismic uncertainties from tomographic velocity inversion and their use in reservoir analysis[END_REF].

Geological realizations obtained by such perturbations of a reference structural model often underestimate structural uncertainties, since interpretation uncertainties are not taken into account. Interpretation uncertainties arise where reflectors cannot be tracked deterministically; this can be a result of an inaccurate imaging velocity model [START_REF] Li | Assessing seismic uncertainty via geostatistical velocitymodel perturbation and image registration: An application to subsalt imaging[END_REF], poor illumination and/or poor resolution [START_REF] Lecomte | 2(3)D convolution modeling of complex geological targets -Beyond 1D convolution[END_REF]. In these situations, multiple geologically possible structural models can be interpreted from the same seismic image [START_REF] Bond | What do you think this is? Conceptual uncertainty in geoscience interpretation[END_REF][START_REF] Bond | Uncertainty in structural interpretation: Lessons to be learnt[END_REF][START_REF] Alcalde | Impact of seismic image quality on fault interpretation uncertainty[END_REF].

Generating multiple structural models from a single data set was one of the earliest methods proposed to evaluate structural uncertainties in geomodeling [START_REF] Thore | [END_REF]. This method opened ways to new challenges, two of which are :

1. Sampling structural uncertainties. This is essentially a structural modeling problem [START_REF] Abrahamsen | Bayesian kriging for seismic depth conversion of a multi-layer reservoir[END_REF][START_REF] Lecour | Modelling of stochastic faults and fault networks in a structural uncertainty study[END_REF]Holden et al., 2003;[START_REF] Wellmann | Towards incorporating uncertainty of structural data in 3D geological inversion[END_REF][START_REF] Wellmann | Validating 3-D structural models with geological knowledge for improved uncertainty evaluations[END_REF][START_REF] Cherpeau | Stochastic structural modelling in sparse data situations[END_REF][START_REF] Julio | Sampling uncertainty about segmented normal fault interpretation using a stochastic downscaling method[END_REF]. In particular, how can we generate multiple structural models respecting both uncertainty bounds and geological rules efficiently? 2. Appraising structural models. This usually amounts to evaluating the model's response for some observed physical phenomenon. In most studies on structural model appraisal, the physical phenomenon that has been used to evaluate the model is fluid flow [START_REF] Suzuki | Dynamic data integration for structural modeling: Model screening approach using a distance-based model parameterization[END_REF][START_REF] Cherpeau | Method for Stochastic Inverse Modeling of Fault Geometry and Connectivity Using Flow Data[END_REF]. Potential fields (gravity and magnetic) have also been used to validate model scenarios [START_REF] Foss | Geologically constrained seismic imaging -Workflow[END_REF][START_REF] Fullagar | Constrained inversion of geologic surfaces -pushing the boundaries[END_REF].

With standard geomodeling workflows, we are able to generate multiple geological scenarios from a seismic image, but we rarely go back to check if these models are consistent with the initial data; doing so would be a way to reduce structural uncertainty. Here, initial data refers to the recorded pre-migrated seismic traces. An interesting approach proposed by [START_REF] Lecomte | Simulated Prestack Local Imaging: a robust and efficient interpretation tool to control illumination, resolution, and time-lapse properties of reservoirs[END_REF] to address this issue is to generate synthetic seismic images from different geological scenarios and compare them to the initial seismic image; a similar approach was adopted by [START_REF] Lallier | Relevance of the stochastic stratigraphic well correlation approach for the study of complex carbonate settings: application to the Malampaya buildup (Offshore Palawan, Philippines)[END_REF]. A major challenge of this approach is that it is not clear how to objectively compare seismic images. An alternative way for reducing interpretation uncertainties in structural modeling is to use "data-driven" interpretation approaches. Here, data-driven interpretation includes all the methods that have been proposed to automatically track horizons and faults, and to build a relative-geologictime function from a seismic image [START_REF] Stark | Relative geologic time (age) volumes-Relating every seismic sample to a geologically reasonable horizon[END_REF][START_REF] Pauget | A global approach in seismic interpretation based on cost function minimization[END_REF][START_REF] Wu | Horizon volumes with interpreted constraints[END_REF][START_REF] Wu | Building 3D subsurface models conforming to seismic structural and stratigraphic features[END_REF]. Data-driven techniques, however, require good quality seismic data [START_REF] Hoyes | A review of global interpretation methods for automated 3d horizon picking[END_REF]; where the geology is complex or lacks reflectivity contrasts, manual interpretation is still necessary.

In this paper, we aim to reduce interpretation uncertainties by appraising structural interpretations using seismic data. Our approach is conceptually similar to that proposed by [START_REF] Lecomte | Simulated Prestack Local Imaging: a robust and efficient interpretation tool to control illumination, resolution, and time-lapse properties of reservoirs[END_REF]: we generate synthetic data from a set of candidate geological models, then we compare the synthetic data against observed data. The candidate geological models are obtained from different interpretations of the same seismic image. There are two main differences between the method proposed here and the method proposed by [START_REF] Lecomte | Simulated Prestack Local Imaging: a robust and efficient interpretation tool to control illumination, resolution, and time-lapse properties of reservoirs[END_REF]. First, we compare initial pre-migrated data instead of seismic images; second, we propose a quantitative comparison instead of a qualitative comparison. The objective of this paper is to investigate how exactly we should compare synthetic data against observed data. We mainly focus on theoretical aspects of the problem and defer the more practical issues for further investigations.

In what follows, we will introduce the problem of appraising structural interpretations using seismic data, and propose a mathematical formulation to describe the process. We then propose some strategies to build a velocity model with structural discontinuities from a structural interpretation and a (smooth) migration velocity model; the resulting velocity model with structural discontinuities is then used to generate synthetic data for the given structural interpretation. Finally, we propose a workflow to design data misfit functions in order to evaluate and rank the different structural interpretations. The proposed method is then applied on a realistic synthetic vertical seismic profiling (VSP) case.

THE STRUCTURAL INTERPRETATION APPRAISAL PROBLEM

We use the model in Figure 1 to illustrate an example of interpretation uncertainty. The model was built by assigning constant velocity values in each major layer of the sandbox model of [START_REF] Colletta | Computerized Xray tomography analysis of sandbox models: Examples of thin-skinned thrust systems[END_REF]. Additional thin-layering structures were added inside each major layer in order to generate more realistic synthetic seismic data [START_REF] Landa | Realistic Finite Differences Modeling -A Case Study[END_REF]. The resulting velocity model was then used to simulate "observed" seismic data using an acoustic, constant density, staggered finite-difference scheme [START_REF] Virieux | SH-wave propagation in heterogeneous media: Velocity-stress finitedifference method[END_REF] with perfectly matched absorbing boundaries [START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]. All subsequent seismic modeling mentioned in this paper were performed with the same code. The reference velocity model in Figure 1 was also used to obtain the migration velocity model in Figure 2a by Gaussian smoothing of the slowness; this migration velocity model was then used to migrate the "observed" data to obtain the seismic image in Figure 2b using Kirchhoff depth migration (e.g. [START_REF] Etgen | An overview of depth migration in exploration geophysics[END_REF].

Figures 3a and3b show two different structural interpretations of the seismic image in Figure 2b from two different interpreters. Both interpretations are geologically possible, but they are different in some regions. Assuming that one is more accurate than the other in the regions where they differ, can we use seismic data to determine which of the two models is more accurate in those regions? To answer this question, we rely on macrolayered velocity models. Macro-layered velocity models are defined as velocity models with discontinuities from structural models and velocity values from the migration model (e.g. Figure 3c and Figure 3d). The migration velocity model is assumed to be kinematically accurate throughout this paper. We use these macro-layered velocity models to compute synthetic data, which are compared against the observed data to appraise the different geological scenarios.

Structural interpretations are typically segments or curves picked by interpreters along structural discontinuities on a seismic image. A structural model is a set of consistent structural surfaces, such as horizons and faults, that represent a geological model [START_REF] Caumon | Surface-Based 3D Modeling of Geological Structures[END_REF]) (e.g. Figure 3a and Figure 3b). In this paper, structural models are built from structural interpretations and there is a one-to-one relation from structural interpretations to structural models. Therefore, appraising structural interpretations is synonymous with appraising structural models herein. The relation from structural models to macro-layered velocity models is assumed to be one-to-one as well.

Let us define a model-misfit as a number that quantifies the mismatch between the reference structural model and a candidate structural model. The reference structural model is the "true" structural model; it is unknown in practice. However, for the sake of argument, we assume the reference model to be known for now. Let us also assume that we can determine this model-misfit value for each model in a given set of candidate structural models. Furthermore, we assume that different models will have different model-misfit values, allowing us to rank those models from the best model to the worst model. This ranking of structural models using model-misfit values will be referred to as model-space ranking of structural models. Objects in the model-space are structural interpretations, or alternatively structural models. We also define a data-misfit as a number that quantifies the mismatch between observed data and the synthetic data generated from candidate structural models. Data-misfit values will allow us to rank structural interpretations in the data-space; we refer to this as data-space ranking. The data-space here refers to the pre-migration domain; i.e., objects in the data-space are seismic traces. We assume that the observed data are the wave equation's response of the reference model (true earth) subjected to some acquisition geometry, just like synthetic data are the wave equation's response of a given candidate (synthetic) model. If the ranking in data-space is the same as the ranking in model-space, then we do not need to know the reference model to determine which among our candidate models are more probable than others. For this reasoning to be applicable, the following conditions have to be satisfied:

1. Different candidate structural models should, overall, have different data-misfit values.

2. Data-misfit values should, overall, be rank-correlated to model-space misfit values.

In practice, a pair of different structural models may exhibit similar seismic response at receivers; in this case, the first condition can be satisfied by ignoring one of the models.

The second condition means that if we plot model-misfits and data-misfits on a scatter plot, we should be able to find a monotonic curve that fits the data.

In general, it is straightforward to design data-misfit functions and compute data misfitvalues. However, it is much harder to design data-misfits that satisfy the conditions mentioned above. The objective of this paper is to investigate what it takes to define data-misfit functions that honor the conditions mentioned above as much as possible, and thus run the process with real data, where the reference (true) model truly is unknown.

MATHEMATICAL FORMULATION

Let M i denote a candidate structural model (e.g. Figure 3a and Figure 3b), the index i taking different values for different interpretation scenarios. Each structural model M i can be used, together with the migration velocity model, to build a macro-layered velocity model, also denoted by M i . Examples of macro-layered velocity models will be presented in the next section (see also Figure 3c and Figure 3d).

Let us introduce the model-space misfit function Φ m i,j that measures the difference between two structural models M i and M j , and denote the mismatch between the reference model and a candidate model

M i as Φ m i = Φ m i,ref .
Similarly, we introduce a data-space misfit Φ d i,j that measures the difference between data computed in two macro-layered models M i and M j , and denote the mismatch between data computed in the macro-layered model M i and the observed (reference) data as

Φ d i = Φ d i,ref . Two models M i and M j are said to be data-different if |Φ d i -Φ d j |
> , for some "small" user-defined threshold > 0.

Let M be a set of data-different structural models. If Φ d i and Φ m i are perfectly rankcorrelated for all models in M , then ranking models in M using Φ d i is the same as ranking models in M using Φ m i . In that case, by definition, the structural models in M can be appraised (i.e. ranked) using seismic data. Therefore, given a set of macro-layered velocity models M , we consider the problem of appraising structural models using seismic data to be solvable in M if the following conditions are satisfied:

Condition 1. |Φ d j -Φ d i | > when i = j. Condition 2. ∞ > Φ d j -Φ d i Φ m j -Φ m i > 0 when i = j.
Given a set of random structural models R, we can build a set of data-different structural models M by progressively moving each model M i from R to M , when M i is data-different from all the models already in M . Therefore, condition 1 can be satisfied in practice. We focus on condition 2 hereafter. Let us define Φ d i (Ω) and Φ m i (Ω) respectively as the dataspace and model-space misfit values corresponding to a specific region Ω in the candidate model M i . The reason for localizing misfits in space will be justified later. In that case, two necessary conditions, although not sufficient, for (the localized version of) condition 2 to be satisfied are:

Condition 3. A change of Φ d i (Ω) ⇒ a change of Φ m i (Ω). Condition 4. A change of Φ m i (Ω) ⇒ a change of Φ d i (Ω).
The symbol ⇒ stands for "implies". Furthermore, suppose that Φ m satisfies the following condition:

Condition 5. A change of Φ m i (Ω) ⇔ a change of the structural model M i in Ω.
Then conditions 3 and 4 become:

Condition 6. A change of Φ d i (Ω) ⇒ a change of the structural model M i in Ω. Condition 7. A change of the structural model M i in Ω ⇒ a change of Φ d i (Ω).
At this point, we have transformed the problem of appraising structural models using seismic data to that of satisfying one condition on Φ m (condition 5) and three conditions on Φ d (conditions 1, 6, and 7). The conditions imposed on Φ d can be checked in practice as they do not depend on the (unknown) reference model. Note that condition 2 implies that we assume the relation between the model and data to be weakly nonlinear, i.e., we assume the perturbations in conditions 3 and 4 to be "small". In other words, we assume, to some extent, that all the candidate structural models are "close" to the (unknown) reference model.

Our goal will be to design Φ d such that conditions 6 and 7 are honored as much as possible. We will argue that for VSP data-misfits, condition 7 is relatively easy to satisfy whereas condition 6 is practically impossible to satisfy in general. As for surface seismic are kinematically equivalent to the migration velocity model, itself assumed to be kinematically accurate, both conditions 6 and 7 can be satisfied in theory for stratified (i.e. layered) models; however, it remains challenging to design and compute misfit function satisfying these conditions. These difficulties make the ideal condition 2 too strict in practice; we replace it with the less ambitious condition:

Condition 8. Φ d is statistically rank-correlated with Φ m .
Condition 2 means that the rank correlation coefficient of the variables Φ m and Φ d is 1, while condition 8 means that it is "high enough".

MACRO-LAYERED VELOCITY MODELS

We propose several options for building macro-layered velocity models from a structural interpretation and a migration model. In this section, we use the illustrative model in

Block macro-layered velocity model

The first type of macro-layered model we propose is obtained by averaging migration slowness values (Figure 4b) in each layer of a structural model. Averaging slowness values preserves traveltimes, as opposed to averaging velocities. The velocity of the i th layer is therefore given by

v i = mean 1 M mig (x)
x in layer i .

(1)

This results in a block velocity model, denoted M b hereafter, as illustrated in Figure 5a.

Figure 5b shows velocity model residuals obtained by subtracting the reference model from the block macro-layered model. Well-log information, when available, can also be used to constrain velocity values assigned in each layer of block macro-layered models. Block macro-layered velocity models are not in general kinematically equivalent to the migration velocity model.

Wire macro-layered velocity model

The second type of macro-layered model is the wire velocity model, M w , illustrated in Figure 6. This model is obtained by

M w (x) = M mig (x) + ∆M b (x), (2)
where ∆M b (x) denotes the velocity contrast at point x; it is equal to zero everywhere except at block boundaries where it is equal to the velocity jump between the velocities in adjacent blocks. Wire macro-layered velocity models are kinematically equivalent to the migration velocity model everywhere expect at geological interfaces.

Reflectivity macro-layered velocity model

We also propose the reflectivity velocity model, M r , illustrated in Figure 7 and obtained by

M r (x) = M mig (x) + M b (x) -Mb (x) , (3)

Density macro-layered velocity models

Finally, we propose a variable density macro-layered model, M d , given by the pair (M mig , M * )

where M * is a density model estimated from M b , by Gardner's law for example [START_REF] Gardner | Formation velocity and density-the diagnostic basics for stratigraphic traps[END_REF]. The variable density model allows us to propose different macro-layered velocity models that are kinematically equivalent to the migration model by proposing different impedance models.

DATA-SPACE MISFIT FUNCTIONS

Seismic data can only contain structural information about subsurface regions that have been illuminated by the wavefield recorded at receivers. Therefore, a data-misfit value computed for a given structural interpretation will only be representative of the illuminated region, and not of the entire model. Furthermore, a misfit value is a scalar, so it summarizes all the interpretation errors in the model in just one value. As a consequence, the best information a global misfit can offer is to determine which of two models is better than the other, but not why (i.e., where the errors are coming from). This motivates the need to look for misfit functions that are localized in space, i.e., data-misfit functions that are a function of space. It is important that we know where the errors are coming from in our model because we might be able to update our structural interpretation in that region to lower the misfit. With this in mind, we propose a data-misfit function of the form

Φ d i (Ω) = s,r x∈Ω ω(s, r, x)χ f o (s, r, t), f i (s, r, t) , ( 4
)
where Ω is the region in the model we would like to evaluate, ω is a weighting function depending on the illumination in Ω, f o (s, r) is an observed seismogram from a shot fired at x s = s and recorded at a receiver at x r = r, f i (s, r) is the corresponding synthetic seismogram computed in the candidate macro-layered velocity model M i , and χ is a data residual function computing the mismatch between observed and synthetic data. We will consider one region of interest Ω at a time; this region can have an arbitrary shape, provided it is illuminated by at least one source-receiver raypath. The data-misfit 4 can be written in the more compact form

Φ d i (Ω) = x∈Ω W χ f o (s, r, t), f i (s, r, t) , ( 5
)
where it is now more apparent that the data-misfit function 5 acts in two steps: first, the residual function χ maps the data plane D × D to a reduced data-space D ; second, the weighting operator W maps the reduced data-space D to the image-space I. Here, the image-space refers to the range of any seismic imaging operator; it is usually referred to as the image-domain. The reduced data-space is, by definition, the range of the residual function χ; a concrete example of a reduced data-space will be presented in the following paragraph. Therefore, if we can find an appropriate linear relation L that maps I to D , the weight function ω can be determined from L * : D → I, where L * = W is either the adjoint or the pseudo-inverse of L. In this paper, we limit ourselves to the case where L * is the adjoint and defer the possibility of L * being the pseudo-inverse for further investigations.

In summary, we propose the following steps for using the localized data-misfit function 5 defined above:

1. Find an appropriate residual function χ : D × D → D .

2. Find an appropriate linear map L : I → D .

Example for vertical seismic profiling data

A common choice of χ : D × D → D for VSP data is the L 2 -norm of data windowed around the first-arrival pick [START_REF] Pratt | Seismic waveform inversion in the frequency domain, part 2: Fault delineation in sediments using crosshole data[END_REF]:

χ f o (s, r, t), f i (s, r, t) = t=τ (s,r)-∆t t=τ (s,r)+∆t f o (s, r, t) -f i (s, r, t) 2 = f (s, r) , (6)
where τ (s, r) is the traveltime from s to r, and ∆t is half the time-window size. In this case, an object in the reduced data-space D is a function of the source and receiver coordinates (s, r); therefore, we need to look for L : f (x) → f (s, r). Let us choose L to be the integration of traveltime along a ray from the source to the receiver:

L[m(x)] = x dl(s, r, x)[m(x)] = τ (s, r) ,
where dl(s, r, x) is a ray-segment centered at x along the ray, and m(x) is the slowness model at x. The adjoint L * is readily available (e.g. Christensen, 2010, chapter 4):

L f (x) , f (s, r) D = s,r x dl(s, r, x)f (x) f (s, r) = x f (x) s,r dl(s, r, x)f (s, r) = f (x), L * f (s, r) I ; (7)
< •, • > I denotes an inner product defined in I. By comparing L * as defined in equation 7 to equation 4, we can identify ω(s, r, x) = dl(s, r, x). Because the rays can always be interpolated such that all ray-segments dl(s, r, x) = 1, we conclude that w(s, r, x) = 1 if x is along the ray connecting s and r,

0 otherwise . ( 8
)
This implies that, for VSP, we can also rewrite the data-misfit 4 as

Φ d i (Ω) = s,r ω(s, r, Ω)χ f o (s, r, t), f i (s, r, t) , (9)
where ω(s, r, Ω) is the length of the ray-segment that actually goes through Ω.

The cartoon in Figure 8 shows the geometrical interpretation of using the data-misfit 4 along with equation 6 and equation 8. The data-misfit is projecting data residuals along rays that illuminate the region of interest Ω; this projection (of data residuals into the imagespace) is equivalent to a generalized Radon transform back-projection (e.g. Toft, 1996, chapter 4). Figure 8a shows rays emanating from sources at the free surface to receivers in a well. The green raypaths illuminate the region of interest Ω, while the black raypaths do not; only the green raypaths will actually contribute to the data-misfit 4. However, the data-misfit is still not localized enough: errors are projected along entire rays since every point x along the green rays has ω(s, r, x) = 1, according to equation 8. Therefore, for χ and ω as given by equation 6 and equation 8, respectively, it is impossible to limit the contribution of the data-misfit 4 to points x ∈ Ω in the most general case. This is a violation of condition 6. Figure 8b shows how we can localize the misfit better by adding an additional well in order to illuminate the region of interest from a different direction.

In this case, points x along the green rays have ω(s, r, x) = 2 if x ∈ Ω (i.e. where the two family of raypaths cross each other), otherwise ω(s, r, x) = 1.

Example for surface seismic data

A common choice of χ : D × D → D for surface data is the L p (p ∈ [1, 2]) difference data residuals [START_REF] Engquist | Application of the Wasserstein metric to seismic signals[END_REF], i.e.

χ f o (s, r, t), f i (s, r, t) = |f o (s, r, t) -f i (s, r, t)| p . (10)
Our numerical experiments suggest that this is not the best choice of χ for appraising structural interpretations using reflection data. Consider for example the reference and interpreted models in Figure 4a and Figure 5a, respectively; reflection shot-gathers from each of these models are shown in Figure 9a and Figure 9b. In Figure 5b, we distinguish interpretation errors from picking errors. Interpretation errors typically result from mispositioning geological interfaces due to poor resolution of a seismic image, while picking errors result from mispositioning geological interfaces due to the fact that it is rarely possible to perfectly track manually a reflector on a seismic image. Looking at the model-space residuals (Figure 5b), picking errors are practically negligible, while in the data-space (Figure 9c) picking errors are as strong as interpretation errors. This is problematic because picking errors are inevitable in practice. In Figure 9c, the picking errors identified by the top green arrow are strong because of the high amplitude of early arrivals; the picking errors identified by the bottom green arrow are important because of the velocity error introduced by the interpretation error above them. This second point motivates the use of macrolayered velocity models that are kinematically equivalent to the migration velocity model when appraising structural interpretations using surface data: if a macro-layered velocity model introduces a velocity error, then all reflectors below this error will be mispositioned in the data-space even if they were positioned at the right place in the image-space by the interpreter. This eventually leads to a violation of condition 6 after back-projecting data errors into the image-space. An illustration of these problems is presented in Figure 10 where we compare individual traces. Figure 10a compares a trace from Figure 9a against a trace from Figure 9b; the difference between these two traces is shown in Figure 10b, where it can be noted that picking errors (δp) are as strong as errors due to mispositioning of horizons (δh). Figure 10c shows a similar experiment but using synthetic data computed in the more kinematically accurate reflectivity macro-layered model in Figure 7; in this case, the additional velocity errors are removed as highlighted by the green circle in Figure 10c.

In addition to being very sensitive to picking errors, equation 10 also suffers from cycleskipping: if a horizon's position is progressively shifted away from its true position, there is a limit beyond which a data-misfit computed using equation 10 will no longer depend on the horizon's position [START_REF] Engquist | Application of the Wasserstein metric to seismic signals[END_REF]. This is a violation of condition 7.

For instance, the data-misfit r, t)] between the two traces in Figure 10a would not change significantly if the mispositioned horizon in Figure 5a was shifted further down because the corresponding reflection-pulse no longer overlaps with the reference one. We therefore conclude that the L p (p ∈ [1, 2]) seismogram difference is not the best choice for χ. In principle, phase data residuals can overcome the aforementioned problems. Here, a phase residual refers to any function of time that measures phase-shifts between two signals. Such residual functions include dynamic-warping-based techniques [START_REF] Hale | Dynamic warping of seismic images[END_REF], and optimal-transport-based techniques [START_REF] Engquist | Application of the Wasserstein metric to seismic signals[END_REF]. The dynamic-time-warping (DTW) phase residual δφ(s, r, t) of f o (s, r, t) and f i (s, r, t) is defined as [START_REF] Hale | Dynamic warping of seismic images[END_REF]:

φ d i (s, r) = t χ[f o (s, r, t), f i (s,
δφ(s, r, t) such that t f o s, r, t -f i s, r, t + δφ(s, r, t) 2 is minimized. ( 11
)
Figure 10d shows the data residuals of the two signals in Figure 10a using this choice of χ. We note that such a misfit is more sensitive to horizon mispositioning errors, and less sensitive to picking errors. However, such a misfit can have some undesired effects as highlighted by the red ellipse. When using equation 10 or 11, an object in the reduced data-space D is a function of time as well as source and receiver coordinates, that is D = D; therefore we need to look for L : f (x) → f (s, r, t). For simplicity, let us choose L to be a Kirchhoff-type modeling operator

L m(x) = x α(s, x, r)S t -τ (s, x, r) m(x) = f (s, r, t) , (12)
where f (s, r, t) is the computed trace recorded at receiver r from a shot at s; S(t) is the source wavelet; α(s, x, r) are "appropriate" migration weights as may be found, for example, in [START_REF] Cohen | Three-dimensional Born inversion with an arbitrary reference[END_REF]; τ (s, x, r) is the traveltime for a ray originating from s scattered off

x and recorded at r; m is the reflectivity model. The adjoint L * is readily available:

L f (x) , f (s, r, t) D =D = s,r,t x α(s, x, r)S t -τ (s, x, r) f (x) f (s, r, t) = x f (x) s,r,t α(s, x, r)S t -τ (s, x, r) f (s, r, t) = f (x), L * f (s, r, t) I . (13)
L * as defined in equation 13 differs from the standard Kirchhoff depth migration (KDM) in that it does not have a derivative along the vertical axis. One often defines [START_REF] Santos | Seismic modeling by demigration[END_REF]; we use this latter definition of KDM hereafter.

KDM f (s, r, t) = ∂ ∂z L * f (s, r, t) (Schneider, 1978), or KDM f (s, r, t) = L * ∂ ∂t f (s, r, t) (14)
By comparing L * as defined in equation 13 to equation 4, we can identify

ω(s, r, x) = t α(s, x, r)S t -τ (s, x, r) . (15)
The interpretation of this result is that the localized data-misfit 4 merely amounts to migrating data residuals for surface data. Figure 11a shows the L * projection of L 1 data residuals of waveforms computed in the model in Figure 4a and waveforms computed in the model in Figure 5a; notable contribution of picking errors (green arrows) and velocity errors (red arrow) are observed. Figure 11b shows the L * projection of L 1 data residuals of waveforms computed in the model in Figure 4a and waveforms computed in the model in Figure 7; velocity errors are dramatically reduced since the model in Figure 7 is more kinematically accurate than the model in Figure 5a. Figure 12a shows the KDM of DTW data residuals of waveforms computed in the model in Figure 4a and waveforms computed in the model in Figure 7. We observe that high-frequency noise is present and that the residual map has two poles: positive contribution above the black line, negative contribution below the black line. Both the high-frequency noise and the bipolarity of projected residuals are a footprint of the derivative along the vertical axis in KDM, and can be removed by using the L * projection as shown in Figure 12b.

APPLICATION & DISCUSSION

Appraising structural interpretations quantitatively is a challenging problem. Consider the eight different structural models in Figure 13 interpreted from the seismic image in Figure 2b by different interpreters. We would like to rank those structural models using VSP data from the most likely to the least likely. For each structural model, we built a block macrolayered velocity model, allowing us to compute VSP synthetic data for each structural model. Observed VSP data were generated in the reference model in Figure 1, with sources covering the free surface and receivers positioned in a vertical well at x = 7500 m. First, we tried to rank the structural models using a data-misfit defined as the L 1 -norm of the entire data set, for each model. Figure 14a shows that the data-space ranking obtained this way is not consistent with the ranking expected in the model-space (the model-space misfits were computed as the L 1 -norm of block macro-layered velocity models and the reference velocity model). We then tried to rank the structural models using the method proposed in this paper. Figure 14b shows the results obtained by limiting the model-misfit and the data-misfit to the region of interest, i.e. the black box in Figure 2b, and by relying on equations 4, 6 and 8 to define the data-misfit. We observe that the scatter plot in Figure 14b satisfies condition 8, while the scatter plot in Figure 14a does not. The model-misfit axis in Figures 14a and14b show that model 1 is the best model. The data-misfit in 14b

was able to indentify model 1 as the best model, whereas the data-misfit in Figure 14a identified model 3 as the best model. In fact, the data-misfit in Figure 14b was able to successfully rank all the models except model 7.

The present paper focuses on misfit functions with the aim of determining the possible reasons as to why the straightforward approach implemented in Figure 14a fails. First, we argue that one should use localized data-misfits Φ d (Ω). It follows that conditions 6 and 7 are necessary conditions for being able to rank structural models using Φ d (Ω), if we expect data-space ranking to be consistent with model-space ranking. These conditions are expected to be valid if the available candidate structural models are "close enough" to the (unknown) reference model. Condition 6 mandates that a perturbation of the data-misfit value Φ d i (Ω) in the region Ω of the i th model should only occur if that model is modified in that specific region Ω. Condition 7 mandates that if the i th model is modified in the region Ω, the data misfit Φ d i (Ω) corresponding to that region should change as well.

For VSP data, condition 7 can readily be satisfied by computing Φ d using only sourcereceiver pairs that illuminate Ω (e.g. using equation 8) and windowing data, along the time axis, around the first arrival (e.g. using equation 6). The difficulty then comes from satisfying condition 6, which is impossible in the most general case as any perturbation of the structural model at any point along any raypath through Ω will lead to a perturbation of Φ d , not just points in Ω. This difficulty can be alleviated, in principle, by adding data that illuminate Ω from different directions as illustrated in Figure 8.

For surface seismic data, using L p (p ∈ [1, 2]) data residuals, we show that velocity errors introduced in macro-layered velocity models can lead to a perturbation of Φ d (Ω) for a region Ω below the velocity error even in the absence of any structural interpretation errors in Ω (e.g. Figure 11a), thereby violating condition 6. We also argue that if a reflector in the structural model is shifted away from its true position until its reflection-pulse no longer overlaps with true reflection-pulse along the time axis, shifting the reflector further in the model would not necessarily affect the data misfit, thereby violating condition 7. However, it is possible, in principle, to approximate condition 7 by using more appropriate data residuals like phase-shift residuals. The challenge then becomes how to compute phaseshifts for complex data sets; this is a subject of ongoing investigations.

The proposed method assumes that the region of interest is well illuminated. We can therefore expect the method to perform poorly for more challenging seismic imaging targets, such as subsalt areas for example. Moreover, appraising structural models using surface seismic data requires a kinematically accurate migration velocity model, which is not easy to obtain in practice. Challenging our approach with a realistic migration velocity model obtained by migration velocity analysis is a subject of ongoing investigations.

CONCLUSION

In this paper, we propose a theoretical formulation and some general solutions to the problem of appraising structural interpretations using seismic data. Assuming that the different structural models are close to the unknown true model, we propose a set of conditions imposed on data-space misfit functions needed for a reliable model appraisal. We argue that misfit functions should be able to localize interpretation errors in the image-space. This localization of errors is achieved by back-projection of data-residuals into the image-space.

It follows that, because it is not possible to truly localize errors using VSP data, one cannot predict which interpretations are more probable than others using VSP data in the most general case. However, it is expected that VSP data can always be used to statistically rank structural interpretations if the data illuminates the target from different directions.

As for appraising structural models using surface seismic data, we expect a better localization of errors, compared to VSP, and therefore a better chance for ranking structural interpretations. The challenge when using surface seismic data is to define appropriate data residual functions. We argue that phase-shift residual functions, such as dynamic-warping and optimal-transport residuals, are good candidates; our current work involves finding robust ways to compute such phase-shift residuals for data acquired in complex geological settings.

FIGURES

Figure

  Figure 4a as the reference model, and the model in Figure 4b as the migration model; the migration model was obtained by smoothing the slowness of the reference model.

where

  Mb is a smooth velocity model obtained by smoothing the slowness of M b . The idea behind the reflectivity macro-layered model comes from the Born approximation which approximates the true velocity model by the sum of a smooth background model plus a reflectivity function (Bleistein et al., 2001, chapter 2); the method proposed here attempts to replace the interpretation-dependent smooth background model Mb with the migration velocity model M mig .

3.

  Find the adjoint L * : D → I. 4. Determine the weight function ω from L * = W.

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 12345 Figure 1: The reference velocity model. It is derived from the sandbox model of Colletta et al. (1991).

Figure 6 :

 6 Figure 6: The same candidate interpretation from Figure 5a represented using a wire macrolayered model (M w ), as described by equation 2.

Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :Figure 12 :FigureFigure 14 :

 78910111214 Figure 7: The same candidate interpretation from Figure 5a represented using a reflectivity macro-layered model (M r ), as described by equation 3.

ACKNOWLEDGMENTS

The authors would like to thank Ludovic Métivier and Pierre Thore for fruitful discussions about the subject. This work also greatly benefited from constructive discussions with Thomas Bodin and Yann Capdeville in the frame of the HIWAI Project. We very thankful to Elizabeth L'Heureux, Isabelle Lecomte, Stan Davis, Xinming Wu for the time they dedicated to review this paper. We would also like to thank Christelle Butault, Francois Bonneau, and Paul Angrand, for some seismic interpretations used in this work.

The traveltime maps needed to evaluate Kirchhoff integral operators were computed using the Madagascar software freely available from www.ahay.org. This work was done in the frame of the RING project at Université de Lorraine. The sponsors of the RING-GOCAD Consortium managed by ASGA are hereby acknowledged for their support. We would also like to thank Paradigm for providing the SKUA-GOCAD software used for structural modeling.


Asmaa Amehraye
email: asmaa.amehraye@enst-bretagne.fr

† D Pastor

Speech enhancement and psychoacoustics

Asmaa Amehraye * , Dominique Pastor †

Abstract

The purpose of this presentation is to describe why and how psychoacoustic models are used to design speech enhancement systems capable of improving speech intelligibility in presence of noise (for mobile communication applications, for instance) or capable of denoising sufficiently well noisy speech signals so as to improve the recognition rate of some automatic speech recognizer (for instance, robot monitoring in noisy environment, hand-free applications on board of vehicles, military fastjets and helicopters).

To begin with, standard methods aimed at denoising speech signals are performed in the spectral domain without taking into account the perceptual characteristics of the speech signal to enhance. They succeed in improving the Signal to Noise Ratio (SNR) but return annoying and unpleasant residual noise known as musical noise.

In the last few decades, psychoacoustic models have then attracted a great deal of interest. The objective is to improve the perceptual quality of the enhanced speech signal. The psychoacoustic model is used to control the enhancement process in order to find the best trade-off between noise reduction, residual noise and speech distortion. The masking phenomenon is the main human auditory system property which is used to design perceptually motivated speech enhancement systems.

The masking phenomenon occurs when the maskee is made inaudible by a stronger signal that occurs simultaneously. We model this phenomenon by a masking threshold calculation that depends on the signal to denoise ( [START_REF] Johnston | Transform coding of audio signals using perceptual noise criteria[END_REF]). Any signal below this threshold is inaudible. The main idea is then to mask partially or totally the distortion. We will present and comment several filters based on this principle ( [START_REF] Hu | Incorporating a psychoacoustic model in frequency domain speech enhancement[END_REF], [START_REF] Amehraye | On the Application of Recent Results in Statistical Decision and Estimation Theory to Perceptual Filtering of Noisy Speech Signals[END_REF]).

Psychoacoustic models are also very popular to assess the quality of speech after coding or enhancement. By using objective tests based on these models, we avoid subjective tests that are very costly in time and resources. The most used objective criteria based on psychoacoustic models are the MBSD (Modified Bark Spectral Distortion) and the PESQ (Perceptual Evaluation of Speech Quality) [START_REF] Quackenbush | Objective Measures of Speech Quality[END_REF], [START_REF] Wang | An objective measure for predicting subjective quality of speech coders[END_REF], [START_REF] Yang | Modified bark spectral distortion measure which uses noise masking threshold[END_REF]. We will present how these objective tests are used and commented in relation with standard and perceptually motivated denoising techniques.


Gilles Pitard
email: gilles.pitard@univ-savoie.fr

Gaëtan Le Goïc

Hugues Favreliere

Serge Samper

Simon-Frédéric Desage

Maurice Pillet

Hugues Favrelière

Discrete Modal Decomposition for surface appearance modelling and rendering

Keywords: Discrete Modal Decomposition, Reflectance Transformation Imaging, Surface appearance, BRDF

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

In many industrial fields, such as micro-mechanics, biomedicine, horology, surface quality is measured visually by human observation of reflected images on the surface, viewed under controlled lighting conditions in order to identify potential defects. The sensory evaluation is often hand-operated, quality controllers vary viewing and light incidence angles on the surface to ensure that they have covered a wide range of configurations. In order to limit the risk of undetected appearance defects, specific configurations are sought (Figure 1): black light e↵ect, white light e↵ect and mirror e↵ect. Based on the observation of this current practice, our approach consists in analysing the interaction of surfaces with their light environment thanks to an estimate of the biangular components of the reflectance function [START_REF] Nicodemus | Directional Reflectance and Emissivity of an Opaque Surface[END_REF][START_REF] Nicodemus | Geometrical considerations and nomenclature for Reflectance[END_REF][START_REF] Koenderink | Bidirectional Reflection Distribution Function expressed in terms of surface scattering modes[END_REF] which quantifies the light reflecting o↵ an object. To adress the modelling of reflectance, three approachs di↵er in complexity and accuracy to approximate the Bidirectional Reflectance Distribution Function (BRDF) which is one of the most complete physical reflection model: L'observateur regarde dans l'axe du rayon spéculaire ( ✓2 ✓1). La surface apparait alors très brillante (cette configuration est aussi appelée e↵et lumière blanche).

L'observateur regarde la surface verticalement et se voit dans le reflet sur la surface (✓2 0 ). • model-based rendering methods consist in the parametrization on analytical BRDF models. This approach is not suitable to adress the appearance quality control due to the complexity of surface textures and light behavior make it di cult to apply to the real materials.

• measuring the complete Bidirectional Reflectance Distribution Function. A dense reflectance sample at a single surface point is generally acquired by a gonioreflectometer. [START_REF] Murray-Coleman | The automated measurement of BRDFs and their application to luminaire modeling[END_REF][START_REF] Dana | Reflectance and Texture of Real-World Surfaces[END_REF][START_REF] Ged | Recognizing real materials from their glossy appearance[END_REF] This method provides a complete description of the interaction light/surface, but will lead to the generation of huge amounts of data and the cost of acquiring data make it incompatible with an application in the context of surface inspection, especially when applied to finding defects.

• image-based rendering methods rely on the acquisition of photometric data (only the direction of illumination varies between images), similar approach as known as Reflectance Transformation Imaging techniques. [START_REF] Macdonald | Visualising an Egyptian Artefact in 3D: Comparing RTI with Laser Scanning[END_REF][START_REF] Zhang | E cient robust image interpolation and surface properties using polynomial texture mapping[END_REF] Thus, our approach is based on the RTI techniques which can accurately model a surface appearance without needing knowledge of its materials or geometrical properties. Several RTI-based techniques have been developped for better measured approximation to reflectance field that uses a set of images -taken from a camera in a fixed position -of the object surface illuminated from di↵erent light directions.

The RTI approach and a sphere capture system are presented in Section 2. After presenting the usual RTI methods (Section 2.

REFLECTANCE TRANSFORMATION IMAGING

RTI approach

A proportional relationship exists between the gray level G of an image and the luminance L( ! u c ) along the optical axis ! u c even if the estimate is biased in presence of cast shadows or overexposed regions. As similar formulations used in numerous papers about Shape-From-Shading, the gray level of an image can be expressed as follows:

9 G = k ⇥ g with g = L( ! u c ) ( 1
)
where k is the factor of proportionality, and g is the normalized gray level. A camera can consenquently be used to make measurements of the luminance of associated pixel values. The RTI principe can divide in 3 steps:

1. Acquisition: The RTI acquisition device provide for each pixel a set of corresponding luminance values for di↵erent angles of illumination. The set of intensities at each pixel location are captured by a fixed camera at position (✓ v , v ) under di↵erent light directions (✓ i , i ).

Modelling:

The bidirectional reflectance function BRDF is approximated from the discrete luminance values to model a continuous information. The surface reflectance is approximated by an RTI method in which di↵erent approaches are detailed from the Section 2.2. The surface reflectance is modelled and represented in the UV texture coordinate system 10 (Figure 2) where (l u , l v ) are the projections of the normalized light vector onto the local basis (L u , L v ) of a particular pixel at the spatial coordinates (u, v) in the studied surface. The descriptors of HSH and DMD techniques are expressed as functions of angles for azimuth and co-latitude ✓ by the following change of variables:

φ View direction! Light direction! Surface Σ! dΣ 1 pixel! (θ v ,φ v ) (θ i ,φ i ) L u L v G L u L v
= arctan 2(l u , l v ) and ✓ = arccos( p 1 l 2 u l 2 v ) (2)

Image acquisition

Figure 3 shows the steps involved in implementing the RTI techniques which require the acquisition of multiple images of a surface with a fixed camera under varying lighting conditions. A system with multiple cameras positioned around a dark sphere was developed. It enables the acquisition of surfaces-of-revolution or freeform surfaces which requires several viewing angles, such as spherical roller bearings (aeronautics), femoral heads of hip implant (medical) or rings (jewellery) because visual coverage of the entire surface can be an issue using a single scene description. The multi-camera system also opens up the creation enriching new methods for 3D reconstruction, surface analysis based on curvatures and slopes, or even measuring transmittance properties of transparent materials. In the present paper, we focus on the upper dome of our spherical structure using only a single camera (cam1) for gathering and processing reflectance transformation imaging. The object is positioned on a glass at the center of the sphere (65 cm diameter). 56 LED light sources are homogenously distributed over the upper hemisphere and above the inspected surface.

Existing modelling techniques 2.3.1 Polynomial Texture Mapping

Polynomial Texture Mappings was originally developped by Tom Malzbender and Al. [START_REF] Malzbender | Polynomial texture maps[END_REF][START_REF] Malzbender | Apparatus for and method of enhancing shape perception with parametric texture maps[END_REF][START_REF] Wolters | Apparatus for and method of converting height fields into parametric texture maps[END_REF][START_REF] Malzbender | Direction-dependent texture maps in a graphics system[END_REF][START_REF] Ritter | System and method for synthesis of parametric texture map textures[END_REF][START_REF] Malzbender | Apparatus for and method of rendering 3D objects with parametric texture maps[END_REF][START_REF] Horton | System and method for combining parametric texture maps[END_REF] The classical PTM model consists in approximating the BRDF from the discrete values of luminance obtained by the device, via polynomial regression. For each pixel, we can seek a polynomial regression 6-vectors of coe cients (a 0 a 5 ) of the quadratic surface defined by the second degree polynomial equation 3. The regressing is carried out using Least Squares formulation (Eq. 4) and solved using Singular Value Decomposition (SVD): Hence, we can interpolate and relight in order to simulate the rendering of a surface under any chosen light direction. PTM descriptors are basic shapes and they provide limited estimation about the overall shape of reflectance field (ie. high-frequency components removed). PTM tends to smooth out the fine details of the surface where sharp shadow boundaries and highly specular surfaces occur.

L(l u , l v ) = a 0 + a 1 l u + a 2 l v + a 3 l u l v + a 4 l 2 u + a 5 l 2 v (3) 2 6 6 6 6 6 6 4 l 2 u0 l 2 v0 l u0 l v0 l u0 l v0 1 l 2 u1 l 2 v1 l u1 l v1 l u1 l

Hemispherical Harmonics

HemiSpherical Harmonics [START_REF] Westin | Predicting Reflectance Functions from Complex Surfaces[END_REF][START_REF] Gautron | A Novel Hemispherical Basis for Accurate and E cient Rendering[END_REF] is an evolution of the PTM technique. This method consists in using a more appropriate projection into a set of basis functions, based on hemispherical harmonics whose shapes are close to reflectance field. HSH have been applied in several topics in computer graphics such as the modelling of BRDFs, irradiance environment maps or invariant object recognition. The hemispherical harmonics H m l are derived from Spherical Harmonic functions (SH) using shifted associated Legendre Polynomials (Eq. 5) to map these functions to a basis set defined only over the upper hemisphere (Eq. 6):

P m l (cos ✓) = P m l (2 cos ✓ 1) and ✓ 2 h 0, ⇡ 2 i ( 5
)
where Legendre polynomials P m l can be computed by applying recurrence relations. [START_REF] Gautron | A Novel Hemispherical Basis for Accurate and E cient Rendering[END_REF] H m l (✓, ) =

:

p 2 Km l cos(m ) P m l (cos ✓) if m > 0 p 2 Km l sin( m ) P m l (cos ✓) if m < 0 K0 l P 0 l (cos ✓) i fm = 0 (6)
where K m l are the hemispherical normalization factors:

K m l = s (2l + 1)(l |m|) 2⇡(l + |m|) (7)
Any surface reflectance function can be decomposed into a HSH series [START_REF] Zhang | E cient robust image interpolation and surface properties using polynomial texture mapping[END_REF][START_REF] Kautz | Fast, arbitrary BRDF shading for low-frequency lighting using spherical harmonics[END_REF][START_REF] Mudge | Image-based empirical information acquisition, scientific reliability, and long-term digital preservation for the natural sciences and cultural heritage[END_REF][START_REF] Krivanek | E↵ects of global illumination approximations on material appearance[END_REF][START_REF] Palma | Dynamic shading enhancement for reflectance transformation imaging[END_REF] of varying order l and degree m and the decomposition residual: where the coe cients C m l can be obtained as the projection of f onto each basis function (Eq. 9):

f (✓v, v ) (✓ i , i ) = n X l=0 l X m= l C m l (✓ v , v )H m l (✓ i , i ) + R n ( 8
)
H 0 0 H 1 -1 H 1 0 H 1 1 H 2 -2 H 2 -1 H 2 0 H 2 1 H 2 2 b)! a)! HSH# DMD# Q 1 Q 2 Q 3 Q 4 Q 5 Q 6 Q 7 Q 8 Q 9 a)! b)!
C m l (✓ v , v ) = Z 2⇡ 0 Z ⇡ 2 0 f (✓ v , v , ✓ i , i )H m l (✓ i , i ) sin ✓ i d✓ i d i ( 9
)
Visual representations of the first few real hemispherical harmonics are given Figure 4 where blue portions represent regions where |H m l | is positive, and yellow portions represent where it is negative.

DISCRETE MODAL DECOMPOSITION

Resulting from vibratory mechanics, Discrete Modal Decomposition (DMD) models a geometrical element in a geometrical form space. DMD is applied in:

1. geometrical tolerancing [START_REF] Favreliere | Modal Tolerancing : From metrology to specifications[END_REF][START_REF] Samper | Modeling of 2D and 3D Assemblies Taking Into Account Form Errors of Plane Surfaces[END_REF] , to visualize the geometrical errors of surfaces and to specify the acceptable variations -form, position and orientation -of the real geometry around the target, 2. multiscale analysis of topographic surface roughness [START_REF] Le Goic | Multi scale modal decomposition of primary form, waviness and roughness of surfaces[END_REF][START_REF] Grandjean | Multi-scalar analysis of hip implant components using modal decomposition[END_REF] , to characterize the form, waviness and roughness defects of a surface, 3. heat source evaluation [START_REF] Pottier | A new filtering approach dedicated to heat sources computation from thermal field measurements[END_REF][START_REF] Pottier | Proposition of a modal filtering method to enhance heat source computation within heterogeneous thermomechanical problems[END_REF] , to estimate the spatial term of a heat di↵usion problem.

Similarly to the Discrete Fourier Transform or the Discrete Cosine Transform, this decomposition enables to make the projection of the measured surface into an eigen basis built from structural dynamics.

Modal basis

The hemispherical domain [0, ⇡ 2 ] ⇥ [0, 2⇡] designates the local hemisphere of BRDFs at surface points where the reflection integral is evaluated. An illuminated ideal di↵use reflecting surface will have equal luminance from all directions which lie in the half-space adjacent to the surface (Lambertian reflectance). The associated geometry (simple half dome in the proposed problem) leads to the definition of the dynamic structural problem which is classically defined as follows:

M.q + K.q = 0 with q = q(✓, )

where M and K stand for the mass and the sti↵ness matrices respectively. Under such formalism, q(✓, ) stands for the displacements which characterizes the modal shapes defined as the harmonic solution of the equation 10.

The eigen vectors Q k with the pulsations w k are solutions of the equation:

q(✓, , t) = +1 X k=1 Q k cos(w k t) (11)
and form the associated modal basis (Q 0 , Q 1 ..., Q n 1 ) where n is the number of modes used for the projection. The point of the present paper is to associate the displacement fields q(✓, , t) to the measured luminances L = (L 0 , L 1 ...., L N ) which draws the reflectance field. The eigen modes are determined by solving the linear system:

(M 1 K 1 w 2 k I)Q k = 0 (12)
PTM (Polynarntal Mappmg) HSH (Hemaphericai Harmonics)

1118=1.11111111 DMD# Q 1 Q 2 Q 3 Q 4 Q 5 Q 6 Q 7 Q 8 Q 9 a)! b)!

Modelling

Thus, any surface reflectance function f can be expressed as the sum of a linear combination of the modal vectors and the decomposition residual R n (Eq.13):

f (✓v, v ) (✓ i , i ) = n X k=1 k (✓ v , v )Q k (✓ i , i ) + R n ( 13
)
where the modal coe cients k can be obtained as the projection of the vector of measured luminances L onto each modes of the non-orthonormal basis (Eq. 14):

k (✓ v , v ) = (Q T k Q k ) 1 Q T k .L (14)
Thus, the modal parametrization uses a basis of adaptative number of modes to better represent surfaces. The first modal vectors composing the modal basis are given Figure 5 with the latter color represention of HSH descriptors. The complexity of shape and frequency increase when the order n increases. If sampled at a high enough rate, the original function can in theory be perfectly reconstructed from the infinite set of samples. In practice, the number of modes used for the projection can not be equal to the number of measured luminances in order to avoid phenomenon known as aliasing or overfitting. So, the number of modes are chosen as to obey the Nyquist-Shannon sampling theorem. A viewer (Figure 6) was implemented in Matlab including PTM, HSH and DMD reconstructions. The modal coe cients can be viewed as a modal amplitude spectrum that approximates the surface points color change as the lighting position changes. We established our RTI format in order to save the modal coe cients encoding the approximative reconstruction of the reflectance spectrum at each pixel in a computer file with the extension .msa (Measuring Surface Appearance). The MSA-viewer allows the controller to move the light in any direction and to aid him in detecting and evaluating particular aspects of surface appearance.

Rendering

Application on high value-added industrial surfaces

The RTI methods are applied to two surfaces of prestige lighter (Figure 7) in which a V-shaped pattern is engraved with an engine turning machine called guilloche to draw vertical lines on a surface at various spacing and widths. The surface can usually be achieved by light hand abrasion after depositing in a thin silver layer.

Aspect defects have been identified on both surfaces, such as strockes and scratches that will be most apparent in the cropped views. After acquiring the sets of 56 images by our system, we compute the three approaches to obtain for each pixel the PTM, HSH and DMD contributions. Taking into account the sampling theorem, we use a modal basis composed of 32 modes, and according to the original embodiment 16 harmonics for HSH and 6 terms for PTM. We can interpolate and relight in the same light direction of the measured images in order to estimate the quality of RTI-based reconstructions whose results are presented in Section 4.2. Our comparative approach focuses on both the space of image representation and the space of lighting because RTI can be regarded as a mapping from the 4D space of lighting and viewing directions to the space of all 2D images. Polynomial Texture Mapping removes the high frequencies of the reflectance field. The smoothing a↵ects the light space, especially in regions where strong specular spikes (high-frequency) are observed (Figure 8). We compute the Peak Signal-to-Noise Ratio (PSNR) which is the most widely used objective image quality metric, in order to assess the quality of reconstruction in both light directions of the device that are chosen arbitrarily. The higher the PSNR, the better the quality of the reconstructed image. Even if the PSNR value indicates the higher quality for DMD reconstruction, it's important to recognize its limitations:

Comparison/discussion

• it relies strictly on numeric comparison based on a pixel-by-pixel comparison of images and does not actually take into account either the location of observable reconstruction degradation nor any level of biological factors of the human vision system, • it does not provide information along the light space of the RTI techniques. PTM technique tends to transform the measured luminances distribution in the di↵use (matte) surfaces (gaussian distribution) while DMD and, to a lesser extent, HSH preserve the shapes of luminance distributions as illustrated in Figure 3.

A consequence of this transformation of luminance distributions can be observed in Figure 9. Our perceptual assessment di↵ers between the image obtained with the light direction 1 and 2. However, changing lighting direction does not much a↵ect the visual appearance of the surface in rendered images produced with the HSH and PTM fitters. DMD can account much luminance variation, for example in low-angled light, the stroke appears as a clear spot in a dark zone like on the original image. Sudden and strong changes in the constrast can automatically draw visual attention to their locations when the user moves the virtual light in any direction.

When the controller has to evaluate the criticity of a detected anomaly, he automatically moves the light source to yield local and global processing of contrast variations.

CONCLUSIONS AND PERSPECTIVES

This paper presents a robust and e↵ective approach to model the local reflectance from the photometric RTI data acquisitions, basing on a parameterization named Discrete Modal Decomposition. By varying the illumination position, the RTI technique aims at enriching the classical information conveyed by images. Each pixel is described by a set of values rather than one value classically; each value corresponding to a specific illumination position. The local reflectance of each pixel can thus be modelled from the set of luminance values associated to each illumination direction. This modelisation is generally performed by applying an interpolation (or approximation) of the discrete acquired data, using a continuous family of functions to describe each pixel local reflectance. Then, this simplified model allows the surface to be virtually rendered for any (chosen) illumination incidence (azimuth and elevation), and the surface can therefore be inspected by varying virtually the direction of illumination. This exploration enhances the micro-geometry and roughness local morphology, and for industrial applications, allows appearance surace quality and perception to be better controlled.

The relevance of the DMD approach is investigated by comparing the results with the ones obtained performing the two most used RTI methods, which are respectively associated to a second order polynomial (PTM) and an hemispherical harmonics (HSH) vector basis. Particularly, this study figures out an interesting feature of the DMD-based reflectance modelisation, that consists in its ability to describe complex local reflectance in terms of angular and intensity luminance variations. If similar results are obtained for lambertzian di↵use reflectances, a more reliable approximation of specular lobes, and more generally of glossy surface local areas is provided by this approach.

1 : 1 -

 11 Angle d'incidence de la source lumineuse 2 : Angle d'observation 1.4.2 Étape d'évaluation 1.4.3 Étape de décision 1 : Angle d'incidence de la source lumineuse 2 : Angle d'observation E↵et lumière noire 2-E↵et brillance 3-E↵et miroir L'observateur regarde la surface verticalement, et l'éclairage est rasant, ce qui donne à la surface un aspect globalement noir, et fait ressortir certaines irrégularités de surface avec un fort contraste (✓1 0).

Tableau 1 . 1 - 1 -

 111 Les 3 e↵ets lumière 1.4.2 Etape d'évaluation Dans cet exemple, on se réfère à un standard écrit pour évaluer et caractériser les anomalies. Les anomalies visuelles sont classées selon 4 grandes catégories, appelées attributs : il s'agit soit d'une marque, soit d'une hétérogénéité, soit d'une déformation ou enfin d'une pollution. L'opérateur doit alors décrire de façon factuelle les anomalies qu'il détecte, en s'appuyant sur le standard défini pour le type de pièce qu'il inspecte. Cet exemple propose d'utiliser 2 types de critères d'évaluation [Baudet 12] : (i) les critères liés aux conditions d'observation. Tableau 1.1 -Les 3 e↵ets lumières 1 : Angle d'incidence de la source lumineuse 2 : Angle d'observation 1.4.2 Étape d'évaluation 1.4.3 Étape de décision 1 : Angle d'incidence de la source lumineuse 2 : Angle d'observation E↵et lumière noire 2-E↵et brillance 3-E↵et miroir L'observateur regarde la surface verticalement, et l'éclairage est rasant, ce qui donne à la surface un aspect globalement noir, et fait ressortir certaines irrégularités de surface avec un fort contraste (✓1 0). L'observateur regarde dans l'axe du rayon spéculaire ( ✓2 ✓1). La surface apparait alors très brillante (cette configuration est aussi appelée e↵et lumière blanche). L'observateur regarde la surface verticalement et se voit dans le reflet sur la surface (✓2 0 ). Tableau 1.1 -Les 3 e↵ets lumière 1.4.2 Etape d'évaluation Dans cet exemple, on se réfère à un standard écrit pour évaluer et caractériser les anomalies. Les anomalies visuelles sont classées selon 4 grandes catégories, appelées attributs : il s'agit soit d'une marque, soit d'une hétérogénéité, soit d'une déformation ou enfin d'une pollution. L'opérateur doit alors décrire de façon factuelle les anomalies qu'il détecte, en s'appuyant sur le standard défini pour le type de pièce qu'il inspecte. Cet exemple propose d'utiliser 2 types de critères d'évaluation [Baudet 12] : (i) les critères liés aux conditions d'observation.Tableau 1.1 -Les 3 e↵ets lumières 1 : Angle d'incidence de la source lumineuse 2 : Angle d'observation 1.4.2 Étape d'évaluation 1.4.3 Étape de décision 1 : Angle d'incidence de la source lumineuse 2 : Angle d'observation

L

  'observateur regarde la surface verticalement et se voit dans le reflet sur la surface (✓2 0 ). Tableau 1.1 -Les 3 e↵ets lumière 1.4.2 Etape d'évaluation Dans cet exemple, on se réfère à un standard écrit pour évaluer et caractériser les anomalies. Les anomalies visuelles sont classées selon 4 grandes catégories, appelées attributs : il s'agit soit d'une marque, soit d'une hétérogénéité, soit d'une déformation ou enfin d'une pollution. L'opérateur doit alors décrire de façon factuelle les anomalies qu'il détecte, en s'appuyant sur le standard défini pour le type de pièce qu'il inspecte. Cet exemple propose d'utiliser 2 types de critères d'évaluation [Baudet 12] : (i) les critères liés aux conditions d'observation.

Figure 1 .

 1 Figure 1. Sensory e↵ects obtained in di↵erent light configurations of a polished lighter (bottom view) by quality controllers

  3), a new modelling technique of angular components of local reflectance from the set of images based on the Discrete Modale Decomposition is presented in Section 3. A comparison with Polynomial Texutre Mappings (PTM) and Hemispherical Harmonics (HSH) is applied to high value-added industrial surfaces in Section 4. It appears that the DMD technique allows us to aid in quantifying particular aspects of surface appearance by preserving spatially local luminance variations in the reconstructed images.

Figure 2 .

 2 Figure 2. Light direction parametrization lu, lv 3. Rendering: The model enables to compute a reflectance function that can be used to interactively display images under novel lighting conditions varying a virtual source light controlled by a user to reveal surface phenomena.

Figure 3 .

 3 Figure 3. Diagram of the automatic acquisition of N images for gathering and processing Reflectance Transformation Imaging

Figure 4 .

 4 Figure 4. Function plots H m l of the Hemispherical Harmonics basis for order l from 0 to 2

Figure 5 .

 5 Figure 5. Function plots Q k of the Discrete Modal Decomposition basis (first 9 modes)

Figure 6 .

 6 Figure 6. MSA-viewer interface. The user can control interactively the light direction by moving a mouse on an intuitive graphic in order to detect and evaluate particular aspects of surface appearance

Figure 7 .

 7 Figure 7. Two bodies of prestige ligther with guilloched pattern

PSNR = 16 Figure 8 .

 168 Figure 8. Cropped views of the original and reconstructed images of lighter 1 and their PSNR value in dB

  -.a+=i :w; !.' / ;.aiil:iL` *aì: = V'T='.-AaL.ilYt.ie .._ '.; f ..Y .!.f--t*4

PSNR = 16 Figure 9 .

 169 Figure 9. Cropped views of the original and reconstructed images of lighter 2 and their PSNR value in dB

Proc. of SPIE Vol. 9525 952523-2 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 11/30/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

Proc. of SPIE Vol. 9525 952523-4 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 11/30/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

Proc. of SPIE Vol. 9525 952523-5 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 11/30/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

Proc. of SPIE Vol. 9525 952523-6 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 11/30/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

Proc. of SPIE Vol. 9525 952523-7 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 11/30/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

Proc. of SPIE Vol. 9525 952523-8 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 11/30/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

Proc. of SPIE Vol. 9525 952523-9 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 11/30/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

Proc. of SPIE Vol. 9525 952523-10 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 11/30/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

Future work will include an assessment of RTI reconstruction quality by analyzing more precisely how this new technique behaves, and particularly, the evolution of luminance local distribution from the acquisition to the RTI modelisation/reconstruction will be investigated. The evaluation of the RTI methods will be carry out statistically in order to better quantify their relevance.


T Hoang
email: tien.hoang@enpc.fr

D Duhamel

G Foret

Dynamical response of a Timoshenko beams on periodical nonlinear supports subjected to moving forces

Keywords: Periodically supported beam, System equivalence, Non-linear, Dynamics, Fourier transform, Timoshenko beam

The periodically supported Timoshenko beam subjected to moving forces has been investigated by numerous researches. The existed models have been developed for linear supports, and this article presents a new one for nonlinear supports. By using a periodic condition and the Fourier series development, the dynamic equation of the Timoshenko beam leads to a relation between the beam displacements and the reaction forces of the supports. This relation does not depend on the support behavior and it exists also for the Euler-Bernoulli beam. Then, the responses can be obtained by combining this relation and the constitutive law of the supports. A numerical method based on discretization of the time and frequency responses has been developed for nonlinear supports.

Moreover, the influence of the beam model has been studied with numerical examples of linear and nonlinear supports. The results show that the Timoshenko beam should be used for the moving forces with high speed and/or the supports with large stiffness.

Introduction

The periodically supported beam subjected to moving forces has been investigated in numerous publications [START_REF] Mead | Free wave propagation in periodically supported, infinite beams[END_REF][START_REF] Mead | Wave propagation in continuous periodic structures: research contributions from southampton[END_REF][START_REF] Metrikine | Vibration of a periodically supported beam on an elastic half-space[END_REF][START_REF] Vostroukhov | Periodically supported beam on a viscoelastic layer as a model for dynamic analysis of a high-speed railway track[END_REF][START_REF] Belotserkovskiy | On the oscillation of infinite periodic beams subjected to a moving concentrated force[END_REF][START_REF] Sheng | Response of infinite periodic structure to moving or stationary harmonic loads[END_REF][START_REF] Sheng | Vibration analysis of the continuous beam subjected to a moving mass[END_REF][START_REF] Nordborg | Vertical rail vibrations: Pointforce excitation[END_REF][START_REF] Nordborg | Vertical rail vibrations: Parametric excitation[END_REF][START_REF] Hamet | Railway noise: Use of the timoshenko model in rail vibration studies[END_REF]. In these articles, the Euler-Bernoulli or Timoshenko beams resting on identical supports at periodical intervals have been considered in steady-state. The response to moving forces are calculated analytically when the supports are linear. However, these models cannot be extended easily for nonlinear supports. Recently, the dynamics of a periodically supported beam has been represented by the system equivalence [START_REF] Hoang | Dynamical response of railway tracks in tunnel[END_REF] by using a periodic condition of reaction forces. This model could work for nonlinear supports, but the author has not presented a method to compute the dynamical responses. Some other researches have considered the model of beams on nonlinear foundations (i.e. the beam is supported continuously) by using the perturbation technique [START_REF] Ansari | Frequency analysis of finite beams on nonlinear KelvinVoight foundation under moving loads[END_REF][START_REF] Kargarnovin | Response of beams on nonlinear viscoelastic foundations to harmonic moving loads[END_REF][START_REF] Ding | Dynamic response of an infinite Timoshenko beam on a nonlinear viscoelastic foundation to a moving load[END_REF], the Galerkin method [START_REF] Ding | Convergence of Galerkin truncation for dynamic response of finite beams on nonlinear foundations under a moving load[END_REF] or the numerical methods [START_REF] Chen | Steady state and stability of a beam on a damped tensionless foundation under a moving load[END_REF][START_REF] Sapountzakis | Nonlinear response of shear deformable beams on tensionless nonlinear viscoelastic foundation under moving loads[END_REF]. This article presents a complete analytic model for the dynamics of beams resting on periodic nonlinear supports. A relation between the beam displacement and the reactions forces has been established from the periodic condition and the dynamic equation of the Timoshenko beam. Then, a numerical method has been developed to compute the response form this relation and the constitutive law of the nonlinear supports. Moreover, a comparison between the Timoshenko and Euler-Bernoulli beam models has been performed with numerical examples of linear and nonlinear behaviors.

Periodically supported Timoshenko beam

Dynamical equations in steady-state

Let's consider an infinite Timoshenko beam resting on identical supports at periodical intervals as shown in Figure 1. The beam is subjected to the moving forces Q j characterized by the distance to the first force D j . Let R n (t) be the reaction force of a support at the coordinate x = nl (with n ∈ Z).

2

R, w ks, ηs f (w)

D j Q 1 Q j v l
Figure 1: Forces applied on a periodically supported beam

In steady-state, we suppose that all supports are equivalent and their responses are described by the same function, but with a delay which equals to the time for the forces to cover the distance between them. In other words, the reaction force repeats when the moving forces pass from one support to another

R n (t) = R t - nl v (1)
where R(t) is the reaction force of the support at the origin of the reference system x = 0. The total force applied on the beam can be represented with the help of the Dirac functions

F (x, t) = ∞ n=-∞ R t - x v δ(x -nl) - K j=1 Q j δ(x + D j -vt) (2)
In addition, we have the dynamic equations of the Timoshenko beam

       ρS ∂ 2 w r ∂t 2 = κSG ∂ 2 w r ∂x 2 - ∂φ r ∂x + F (x, t) ρI ∂ 2 φ r ∂t 2 = EI ∂ 2 φ r ∂x 2 + κSG ∂w r ∂x -φ r (3)
where w r , φ r are the displacement and the rotation of the section of the beam;

ρ, E are the density, the Young's modulus and S, I, κ, G are the section, the 30 inertia, the shear coefficient and the shear modulus of the beam.

Equations ( 2) and (3) define the dynamics of the beam and its supports in the steady-state. Thereafter, we will resolve these equations by performing a Fourier transform with regard to the time t, and then using the Fourier series development with regard to x. Let's denote ∂ t , ∂ x the derivations with regard to t and x. By performing the Fourier transform of equation ( 3) with regard to t, we obtain

     κSG ∂ x φr = κSG ∂ 2 x ŵr + ρSω 2 ŵr + F -κSG ∂ x ŵr = EI ∂ 2 x φr -(κSG -ρIω 2 ) φr (4)
where the hat stands for the Fourier transform with regard to t. Particularly, we obtain the following result from equation ( 2)

F = e -iω v x   R(ω) ∞ n=-∞ δ(x -nl) - K j=1 Q j v e -iω v Dj   (5)
Thus, e iω v x F is periodic with regard to x. Therefore, if we put

ŵr = Ψ(x, ω)e -iω v x
and φr = Φ(x, ω)e

-iω v x (6)
equation ( 4) becomes

     κSG ∂ x Φ -iω v Φ = κSG ∂ 2 x Ψ -2 iω v ∂ x Ψ -ω 2 v 2 Ψ + ρSω 2 Ψ + e iω v x F κSG iω v Ψ -∂ x Ψ = EI ∂ 2 x Φ -2 iω v ∂ x Φ -ω 2 v 2 Φ -(κSG -ρIω 2 )Φ (7)
By the Floquet's theorem [START_REF] Walter | Ordinary differential equations[END_REF], equation [START_REF] Sheng | Vibration analysis of the continuous beam subjected to a moving mass[END_REF] has a periodic solution. We can find this solution by using the Fourier series developments of Φ and Ψ (see Appendix A). Thereafter, by combining the results of Ψ and Φ with equation ( 6) we obtain

         ŵr (x, ω) = R(ω) ∞ n=-∞ pn e ix( 2πn l -ω v ) -p0 l v K j=1 Q j e -iω v (x+Dj ) φr (x, ω) = R(ω) ∞ n=-∞ qn e ix( 2πn l -ω v ) -q0 l v K j=1 Q j e -iω v (x+Dj ) (8)
where pn , qn (n ∈ Z) are the Fourier coefficients of Ψ, Φ calculated by equation (A7) in Appendix A. We can reduce equation ( 8) by defining η(x, ω), γ(x, ω) as

follows          η(x, ω) = ∞ n=-∞ pn e ix( 2πn l -ω v ) γ(x, ω) = ∞ n=-∞ qn e ix( 2πn l -ω v ) (9)
Indeed, η, γ in the last equation can be reduced to simple analytical functions as shown in equations (B11) and (B14) of Appendix B. Then, by substituting equation ( 9) into equation [START_REF] Nordborg | Vertical rail vibrations: Pointforce excitation[END_REF], we obtain

         ŵr (x, ω) = R(ω)η(x, ω) -p0 l v K j=1 Q j e -iω v (x+Dj ) φr (x, ω) = R(ω)γ(x, ω) -q0 l v K j=1 Q j e -iω v (x+Dj ) (10)
Equation ( 10) is a simple relation between the Fourier transforms of the beam displacement and the reaction force. This is a result of the periodicity condition and the dynamic equation of the Timoshenko beam, which do not depend on the support behaviour. Once the reaction force is calculated, this equation can 35 be used to compute the response of the beam. In the next section, we will introduce a system equivalence based on this relation.

System equivalence

In order to calculate the reaction force of the support, we need to compute the displacement of the beam at the support position w(t) = w r (0, t), or its Fourier transform ŵ(ω) = ŵr (0, ω). By substituting x = 0 into equation ( 10),

we have

ŵr (0, ω) = R(ω)η(0, ω) -p0 l v j Q j e -i ω v Dj (11)
Hence, we can also write

R(ω) = K T ŵ(ω) + Q T ( 12
)
where K T and Q T are defined by

K T = η -1 (0, ω) and Q T = K T p0 l v K j=1 Q j e -i ω v Dj (13)
with p0 (ω), η(0, ω) are calculated by equations (A10) and (B15) in the appen-

40 dices.
Equation ( 12) is a linear relation between the force and the displacement applied on the support at x = 0, and it holds for all supports because of the periodicity condition. This relation is the same as the constitutive law of an equivalent spring with stiffness K T and pre-force Q T . Therefore, we call the system equivalence of a periodically supported Timoshenko beam, which existed also for Euler-Bernoulli beams (see [START_REF] Hoang | Dynamical response of railway tracks in tunnel[END_REF]). Equation ( 12) explains the distribution mechanism of the moving forces Q j to the supports via the beam. When a moving force comes toward and leaves away the support along the direction of the beam, the reaction force of the support increases and decreases respectively.

This process is the same as a force applied on the support via the equivalent spring.

The comparison of the system equivalences with two parameters stiffness and pre-force for Euler-Bernoulli and Timoshenko beams is shown in Table 1.

We see that the stiffness K depends on two parameters λ 1,2 and C 1,2 which are different between the two beam models. However, if the shear modulus κG and the ratio E/ρ tend to infinity, λ 1,2 and C 1,2 of the Timoshenko beam tend to the ones of the Euler-Bernoulli beam. Therefore, the stiffness and pre-force for the two beam models are equivalent when the Timoshenko beam does not include the shear modulus κG and the ratio E/ρ. This phenomenon agrees well with the beam theories.

R(ω) = K(ω) ŵ(ω) + Q(ω) Q(ω) = K(ω) p0 l v K j=1 Q j e -iω D j v K(ω) = 2EI(λ 2 1 + λ 2 2 ) C1 λ1 sin lλ1 cos lλ1-cos lω v -C2 λ2 sinh lλ2 cosh lλ2-cos lω v -1
Euler -Bernoulli beam Timoshenko beam In other words, the influence of the beam models is more important at high frequencies. In the next sections, we will calculate the responses for two beam 70 models with linear and nonlinear support behaviours in order to compare the two beam models.

λ 2 1,2 = ρSω 2 EI C 1,2 = 1 p0 l = 1 EI ω 4 v 4 -ρSω 2 λ 2 1,2 = ρSω 2 EI + ω 4 4 ρ E -ρ κG 2 ± ω 2 2 ρ E + ρ κG C 1,2 = 1 - ρIω 2 ∓EIλ 2 1,2 κSG p0 l = κSG-ρIω 2 +EI ω 2 v 2 κSG EI ω 4 v 4 -ρSω 2 -ρSI(κG+E-ρv 2 ) ω 4 v 2

Calculation of responses

Linear supports

Let's consider the viscoelastic supports represented by a stiffness k s and a damping coefficient η s . The dynamical stiffness of the supports is K s = k s +iωη s and their constitutive law is given by By combining equations ( 12) and ( 14), we obtains the response of the supports

R(ω) = -K s ŵ(ω) (14)
       R(ω) = K s Q T K T + K s ŵs (ω) = -Q T K T + K s (15)
Then, by substituting the last results and equation [START_REF] Kargarnovin | Response of beams on nonlinear viscoelastic foundations to harmonic moving loads[END_REF] into equation [START_REF] Hamet | Railway noise: Use of the timoshenko model in rail vibration studies[END_REF] we obtain the response of the beam

       ŵr (x, ω) = Q T K s η(x, ω) K s + K T -η(0, ω)e -i ω v x φr (x, ω) = Q T K s γ(x, ω) K s + K T -γ(0, ω)e -i ω v x (16)
Equations ( 15) and ( 16) describe completely the responses in the frequency The responses have been calculated for the linear support by using equation [START_REF] Chen | Steady state and stability of a beam on a damped tensionless foundation under a moving load[END_REF] with the beam parameters given in table 2 and the support parameters 80 k s = 20 MN/m, η s = 0.1 MNs/m. Figure 4 shows the reaction force and the displacement of the support in the frequency domain and the time domain.

Although the equivalent stiffness of the two beam models are different, the responses are almost the same because they are affected more by the lower frequencies. Moreover, the larger velocity leads to responses of smaller amplitudes, but of greater frequency range. Next, we will investigate the influence of the beam model on the responses by parametric studies.

Influence of beam models

The Timoshenko beam theory shows that the two beam models can be equivalent in static when EI κl 2 SG ≪ 1. It is remarkable that this value corresponding to the previous example is 0.0715, which is closed to the limit of this condition. Now we consider a system of supports that has a nonlinear behaviour as shown in figure 1. The reaction force of a support is given by

R(t) = -(k s w + η s w ′ + f (w)) (17)
Combining equations ( 12) and (17) gives a system of two equations with two unknowns R and w. In order to simplify the linear part of the constitutive law, we put

P (t) = R + k s w + η s w ′ (18)
By performing the Fourier transform, we can also write as follows

P (ω) = R + K s ŵ ( 19
)
where K s = k s + iωη s . By combining equations ( 19) and ( 18) with equations ( 12) and ( 17) respectively, we obtain

     P = -f (w) P = (K T + K s ) ŵ + Q T (20)
When f (w) = 0, the last equation is exactly the same result as for linear supports presented in equation [START_REF] Ding | Convergence of Galerkin truncation for dynamic response of finite beams on nonlinear foundations under a moving load[END_REF]. When f (w) = 0, we got one equation in the frequency domain and the other in the time domain. To resolve this problem, we will discretize the unknown functions P (t) and w(t) and their Fourier transforms. It is remarkable that these functions (and their Fourier Therefore, we can use the Shannon sampling theorem [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] to discretize these functions.

Numerical method

Suppose that it exists t max , ω max provided that

     P (t) ≃ 0, w(t) ≃ 0 if |t| ≥ t max P (ω) ≃ 0, ŵ(ω) ≃ 0 if |ω| ≥ ω max (21)
We discretize the responses w(t), P (t) and their Fourier transforms by considering their values at t = nT and ω = nΩ with |n| ≤ N . The discretization has to satisfy the Shannon sampling theorem [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]

ω max ≥ π T and t max ≥ π Ω (22)
Because t max = N T and ω max = N Ω, the last condition becomes

N T Ω ≥ π (23)
The discretization of the first equation of (20) becomes

P n = -f (w n ) |n| ≤ N (24)
where P n = P (nT ) and w n = w(nT ). In order to discretize the second equation of (20), we need to compute the Fourier transforms of the discretized responses P n and w n , which are given by (see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF])

       Pd (ω) = n P n e -inT ω = P (ω) T |ω| ≤ π T ŵd (ω) = n w n e -inT ω = ŵ(ω) T |ω| ≤ π T (25)
Therefore, the discretization of the Fourier transforms ( P , ŵ) is deduced from the last equation

       Pm = P (mΩ) = T n P n e -inmT Ω |mΩ| ≤ π T ŵm = ŵ(mΩ) = T n w n e -inmT Ω |mΩ| ≤ π T (26)
By combining the conditions of equations ( 23), ( 24) and the last equation, we obtain

N ΩT = π (27)
and we have

       Pm = T N n=-N P n e -inmT Ω |m| ≤ N ŵm = T N n=-N w n e -inmT Ω |m| ≤ N (28)
We can also rewrite the last equation in vector forms P = A P and ŵ = A w

where P , w are column vectors of P n , w n (-N ≤ n ≤ N ) and A is a square matrix defined from equation (28). Then, by substituting the last equation into the second equation of (20), we obtain

A P = K A w + Q (30)
where Q is a column vector of Q n = Q T (nΩ) and K is a diagonal matrix of 110

K n = K T (nΩ) + K s (nΩ) (with -N ≤ n ≤ N ).
Finally, by combining equations ( 24) and (30), we obtain

K A w + A f + Q = 0 ( 31
)
where f is a column vector of

f n = f (w n ).
Equation ( 31) is nonlinear with regard to the unknown w. If f (w) ≡ 0 (linear support), we can obtain the solution of this equation as follows

w 0 ≡ -K A -1 Q = -A -1 K -1 Q (32)
The last result is exactly the discretization of the analytic solution in equation [START_REF] Ding | Convergence of Galerkin truncation for dynamic response of finite beams on nonlinear foundations under a moving load[END_REF]. In a general case, we can use a numerical solver for equation (31). In this article, we use the solver "fsolve" of MATLAB . A difficulty is that this solver does not compute complex solutions. We need to decompose the unknown w into their real and imaginary parts. The initial value of the solver is taken by the linearised solution presented in equation (32).

Example

Let's consider a nonlinear support described by the constitutive law in equation ( 17 shows the results with the two beam models for linear and nonlinear behaviours.

We see that the responses have similar forms for the two beams models and support behaviours. In particular, the reaction forces in the frequency domain have the same maximum pick at zero, but the nonlinear support has larger amplitudes in the high frequencies.

Conclusion

The system equivalence of a periodically supported beam has been developed

for the Timoshenko beam and applied to nonlinear supports. This model has been compared with the existing model for an Euler-Bernoulli beam. Then, a numerical method has been developed for nonlinear supports based on the time-frequency transform which permits a fast calculation of the response of the support systems. Moreover, the numerical examples showed that the influence of the beam model is more important when the speed of the loads is higher and/or the support stiffness is larger.

Remark : We can obtain the expressions of p0 , q0 by substituting n = 0 into equation (A7)

       p0 = κSG-ρIω 2 +EI ω 2

FrequencyFigure 2 :FrequencyFigure 3 :

 23 Figure 2: Equivalent stiffness of Timoshenko (circle) and Euler-Bernoulli (continuous line) beams with moving velocity 150 km/h (green) and 300 km/h (red)

75Figure 4 :

 4 Figure 4: Responses of the supports with Timoshenko (circle) and Euler-Bernoulli (continuous line) beams with a moving velocity 150 km/h (green) and 300 km/h (red)

  Therefore, it is necessary to study the influence of the beam model on the dynamic responses. The external parameters which could affect on the responses are the moving force speed and the stiffness of the support. Figures5 and 6

Figure 5 :Figure 6 :

 56 Figure 5: Influence of the moving force speed on the reaction forces

  105 transforms) have bounded supports (see the numerical examples in section 3.1).

6 FrequencyFigure 7 :

 67 Figure 7: Responses of the supports with Timoshenko (circle) and Euler-Bernoulli (continuous line) beams with a linear (green) and a nonlinear (red) support behaviours

v 2 l κSG EI ω 4 v 4 -ω 4 v 4 - 4 v 2 ( 1 sin λ 1 lC 2 λ 2 sinh λ 2 l

 244421122 ρSω 2 -ρSI(κG+E-ρv 2 ) ρSω 2 -ρSI(κG+E-ρv 2 ) ω A10) Appendix B: Analytical expressions of η(x, ω) and γ(x, ω) 145 The function z n in equation (A8) is a second-order polynomial with regard to ω v -2πn L 2 and we can rewrite it as follows z n (ω) = κSGEI ω vsimplify the expression of pn , we substitute equation (B1) into the first equation of (A7), we can write pn (taking the equality of the coefficients of two expression of pn in equations (A7) and (B3), we obtain deduce the analytical expression of η(x, ω) in equations (9), we need to calculate the infinite sum of the series which are presented in equation (B5). Let's consider the following functionf (x) = sin λ 1 (l -x) + e -i ωl v sin λ 1 x cos lλ 1 -cos ωl v e i ω v x ∀x ∈ [0, l](B6)Remark : We obtain the expression of η(0, ω) by substituting x = 0 cos lλ 1 -cos ωl v cosh lλ 2 -cos ωl v (B15)

Table 1 :

 1 System equivalence for Euler-Bernoulli and Timoshenko beams

Table 2 :

 2 Parameters of a periodically supported beam

	Mass density	ρ kg/m 3	7850
	Young's modulus	E GPa	210
	Shear modulus	G GPa	81
	Shear coefficient	κ		0.4
	Section	S m 2	7.69e-3
	Flexion inertia	I	m 4	3.055e-5
	Sleeper spacing	l	m	0.6
	Moving force	Q kN	100
	Speed of moving force v	km/h	300

Example

Figures 2 and 3 show an example of the stiffness and the pre-force for the Timoshenko and the Euler-Bernoulli beams with the parameters presented in

Table 2 .

 2 We see that the two beam models give almost the same equivalent pre-

65

force. Otherwise, the Timoshenko beam gives a smaller equivalent stiffness than the Euler-Bernoulli beam. It is remarkable that the difference takes place only at the maximum peaks of the stiffness which correspond to high frequencies.

Acknowledgments

This work has been developed in the context of a partnership between Eurotunnel Group and Ecole des Ponts ParisTech. The authors would like to thank the personnel of Eurotunnel Group for their support, specially J.L. Pochet, P. Joyez and R. Caby.

Appendix A: Calculation of Ψ(x, ω) and Φ(x, ω)

We will find the periodic solution of equation [START_REF] Sheng | Vibration analysis of the continuous beam subjected to a moving mass[END_REF] q n (ω)e i2πn x l (A1) with p n , q n are computed by

It is remarkable that the Fourier coefficients of e i ω v x F in equation ( 5) are calculated by

where δ 0n = 1 if n = 0 and δ 0n = 0 if n = 0. Therefore, by performing the Fourier series development of equation [START_REF] Sheng | Vibration analysis of the continuous beam subjected to a moving mass[END_REF] and by using the last equation, we obtain

Then, the following result is deduced from the second equation of (A4)

Then, by substituting the last equation into equation (A4), we obtain

where pn , qn are calculated by

with

By substituting equation (A6) into equation (A1), we obtain:

Because f (0) = f (l), the Fourier coefficients of this function are calculated by

We can write this function in the Fourier series expansion as follows

By substituting equations (B6), (B7) into the last equation, we obtain

Similarly, we have

Finally, the analytical expression of η(x, ω) in equations ( 9) can be obtained by combining (B3) and the two last results

In the similar way, we can deduce the analytic expression of γ(x, ω) by applying the derivation with regard to x of equations (B9) and (B10)

Then, by combining the two last equations and equations (B4), we obtain the analytical expression of γ(x, ω) in equation ( 9)

cos λ 1 (l -x) -e -i ωl v cos λ 1 x cos lλ 1 -cos ωl v cosh λ 2 (l -x) -e -i ωl v cosh λ 2 x cosh lλ 2 -cos ωl v (B14)