Julien Réthoré
email: julien.rethore@lamcos.insa-lyon.fr

Stéphane Roux
email: stephane.roux@lmt.ens-cachan.fr

François Hild
email: francois.hild@lmt.ens-cachan.fr

Optimal and noise-robust extraction of Fracture Mechanics parameters from kinematic measurements

Keywords: Full field measurement, noise sensitivity, stress intensity factors, uncertainty

The paper is devoted to an optimal (i.e. noise robust) determination of stress intensity factors and crack tip locations based on a displacement field measured over an arbitrarily shaped domain. As the minimization of the noise sensitivity is included within the proposed extraction technique, this is especially dedicated to corrupted displacement fields, e.g. as measured by an optical technique. The main idea is to construct for mode I and II fields an extracting function so that its L 2 scalar product with the actual displacement field yields the sought parameter. The extracting function is also constrained to be orthogonal to a set of admissible elastic fields. Two applications are considered to illustrate the technique. The first example deals with a fatigue crack in steel for which small scale yielding occurs. A second example with a

One of the central issues of Fracture Mechanics is to evaluate global parameters (e.g. stress intensity factors) from a given displacement field. This is an important question for numerical modeling [START_REF] Parks | The virtual crack extension method for nonlinear material behavior[END_REF][START_REF] Attigui | Mixed-mode separation in dynamic fracture mechanics: New path independent integrals[END_REF][START_REF] Goz | Domain integral formulation for stress intensity factor computation along curved three-dimensional interface cracks[END_REF][START_REF] Rajaram | Application of domain integral methods using tetrahedral elements to the determination of stress intensity factors[END_REF][START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets. Part I: Mechanical model[END_REF] and also an essential problem for experimental characterization [START_REF] Mcneill | Estimation of stress intensity factor by digital image correlation[END_REF][START_REF] Huntley | Measurement of crack tip displacement field using laser speckle photography[END_REF][START_REF] Forquin | A Method to Determine the Toughness Scatter of Brittle Materials[END_REF][START_REF] Réthoré | Estimation of mixed-mode stress intensity factors using digital image correlation and an interaction integral[END_REF][START_REF] Roux | Stress intensity factor measurements from digital image correlation: Post-processing and integrated approaches[END_REF][START_REF] Yoneyama | Evaluating mixed-mode stress intensity factors from full-field displacement fields obtained by optical methods[END_REF] of cracks. One of the main difficulties of this problem is the issue of scale, namely, not too large (to avoid structural effects) but not too small (to keep clear from the inelastic behavior, or numerical discretization biases).

An additional difficulty arises when the displacement field is extracted from experimental data that are corrupted by measurement uncertainties. Minimizing the effect of this noise on the determination of Fracture Mechanics parameters is a key property, which is addressed in the present study.

Two approaches of Linear Elastic Fracture Mechanics are usually distinguished, namely, the energy and the kinematic (or static) approaches. The energy release rate G [START_REF] Griffith | The Phenomenon of Rupture and Flow in Solids[END_REF] and J-integral [START_REF] Rice | A Path Independent Integral and Approximate Analysis of Strain Concentrations by Notches and Cracks[END_REF] and are two complementary methods that evaluate the amount of energy released in a virtual extension of the crack. The kinematic approach defines the stress intensity factors [START_REF] Irwin | Analysis of the Stresses and Strains near the End of a Crack Traversing a Plate[END_REF] for a 2D solid as a limit

K I ∝ lim r→0 µ 2π r u 2 K II ∝ lim r→0 µ 2π r u 1 (1)
In 2D linear elasticity, the following result

J = G ∝ K 2 I + K 2 II µ (2)
unifies the two approaches. One of the properties of these quantities is their quadratic dependence in the displacement field and the stress intensity factors respectively. From this observation, a class of methods devoted to mixed mode separation has been developed (see, e.g. Refs. [START_REF] Parks | A stiffness derivative finite element technique for determination of crack tip stress intensity factors[END_REF][START_REF] Chen | Conservation Laws in Elasticity of the J-Integral Type[END_REF][START_REF] Bui | Associated path independent J -integral for separating mixed modes[END_REF]). Among those, the interaction integral reads as a trilinear form of the actual displacement field u, an auxiliary field u a and a virtual crack extension (VCE) field q

I int = - Ω σ a ml u m,l δ kj -(σ a ij u i,k + σ ij u a i,k ) q k,j dΩ (3)
where Einstein's convention is used, The auxiliary field is chosen to extract separately mode I and II stress intensity factors. Usually, the asymptotic crack tip fields in modes I and II are used to obtain K I and K II , respectively.

The VCE field is constrained by a few assumptions. First, it must vanish outside Ω, be equal to x 1 at the crack tip and be parallel to x 1 along the crack faces. For a given choice of q, I int is a bilinear symmetric positive definite form. Thus the interaction integral may be viewed as a scalar product related to the energy norm defined by J or G. The stress intensity factors are then estimated from this energy-based scalar product between the actual and auxiliary displacement fields.

The most popular stress intensity factor extracting methods dedicated to experimental applications is based on least squares minimization of the gap between the asymptotic crack tip field and the actual displacement field. One focuses on the crack opening displacement only [START_REF] Forquin | A Method to Determine the Toughness Scatter of Brittle Materials[END_REF] or on full fields [START_REF] Mcneill | Estimation of stress intensity factor by digital image correlation[END_REF][START_REF] Roux | Stress intensity factor measurements from digital image correlation: Post-processing and integrated approaches[END_REF][START_REF] Yoneyama | Evaluating mixed-mode stress intensity factors from full-field displacement fields obtained by optical methods[END_REF].

The interaction integral has recently been used to measure stress intensity factors [START_REF] Réthoré | Estimation of mixed-mode stress intensity factors using digital image correlation and an interaction integral[END_REF] and optimized [START_REF] Réthoré | Noise-robust Stress Intensity Factor Determination from Kinematic Field Measurements[END_REF] to reduce the influence of measurement uncertainties on the determination of the sought parameter. The reduction of the noise sensitivity is also one of the aims of the approach proposed herein.

In the present paper, it is shown that all these techniques, interaction integrals as well as least squares extraction, can be recast within the framework detailed in Section 2. It will be shown that the optimized version of the interaction integral [START_REF] Réthoré | Noise-robust Stress Intensity Factor Determination from Kinematic Field Measurements[END_REF] and least squares extraction are local optima with respect to noise sensitivity. The technique proposed herein allows for the numerical determination of an optimal extracting function. Theoretical crack displacement fields and corresponding Fracture Mechanics parameters are detailed in Section 2 as well as the formulation of an optimal noise-robust extractor. A numerical resolution is also proposed in this section. Then the Digital Image Correlation (DIC) algorithm used herein to measure displacement fields is briefly presented to derive its noise sensitivity in Section 3. such that z = re iθ . The crack tip is located at the origin, z = 0, and the semi-infinite crack extends in the direction θ = ±π.

Last, in

The displacement field u is conventionally represented by its complex writing, u = u 1 + iu 2 . It is expanded as a double series

u(r, θ) = i=I,II n a n i φ n i (r, θ) (4)
where a n i are the sought parameters, and φ n i reference fields that form the complete set of linear elastic fields satisfying a condition of zero traction along the crack path:

• even n

φ n I (r, θ) = r n/2 κe inθ/2 - n 2 e i(4-n)θ/2 + ( n 2 + 1)e -inθ/2
(5)

φ n II (r, θ) = ir n/2 κe inθ/2 + n 2 e i(4-n)θ/2 -( n 2 -1)e -inθ/2 (6)
• odd n

φ n I (r, θ) = (-1) (n+1)/2 r n/2 κe inθ/2 - n 2 e i(4-n)θ/2 + ( n 2 -1)e -inθ/2 (7)
φ n II (r, θ) = i(-1) (n-1)/2 r n/2 κe inθ/2 + n 2 e i(4-n)θ/2 -( n 2 + 1)e -inθ/2 (8)
where κ is Kolossov's constant, κ = (3 -ν)/(1 + ν) for plane stress or κ = (3 -4ν) for plane strain conditions.

The partition between mode I and II is performed on the basis of the following symmetries (no discontinuity is involved here since 2 is even). Larger indices correspond to "subsingular" fields that may capture the remote heterogeneity of the loading, but do not affect the fracture behavior at the crack tip.

φ n I (z) = φ n I (z) φ n II (z) = -φ n II (z) (9)
Classically, i.e. within the framework of Linear Elastic Fracture Mechanics, only n ≥ 0 fields are considered, since negative indices would correspond to physically unrealistic displacement fields with a diverging elastic energy density at the crack tip. However, inelastic behaviors are typically activated in the vicinity of the crack tip, in the so-called process zone, so that the above linear elastic decomposition does not hold. If a small domain close to the crack tip that contains all non-linearities is omitted, the displacement field in the complementary domain will naturally contain contributions from such n < 0 fields [START_REF] Hamam | Stress intensity factor gauging by digital image correlation: Application in cyclic fatigue[END_REF][START_REF] Réthoré | An extended and integrated digital image correlation technique applied to the analysis fractured samples[END_REF][START_REF] Henninger | Enriched kinematic fields of cracked structures[END_REF] (referred to as "supersingular" in the sequel).

They are necessary to account for the most general problem where both the small and large r regions are omitted from the elastic description.

Fracture Mechanics parameters

Considering the expression of the general crack tip field in Equations [START_REF] Chen | Conservation Laws in Elasticity of the J-Integral Type[END_REF] to [START_REF] Griffith | The Phenomenon of Rupture and Flow in Solids[END_REF], the following recurrence formulas are obtained

∂φ n j ∂x 1 = - n 2 φ n-2 j ( 10
)
The derivative with respect to the crack tip position x c being the opposite of the above equation, it shows that a mis-positioning of the crack tip along its axis will give rise to a supersingular contribution φ -1 j , away from the immediate vicinity of the crack tip. Therefore, in order to estimate the crack tip position, the displacement field is decomposed by using n = 1 and n = -1 order functions

a 1 j φ 1 j + a -1 j φ -1 j = a 1 j φ 1 j + 2 a -1 j a 1 j ∂φ 1 j ∂x c (11)
For the symmetric part of u, i.e. for mode I, this result is interpreted as a

Taylor expansion of the usual crack tip field whose tip would be shifted from the actual one by an offset distance x I such that

x I = 2 a -1 I a 1 I (12)
Calculating x I from a -1 I (i.e. ∝ K -1 I ) and a 1 I (i.e. ∝ K I ), an estimate of the offset between the assumed crack tip position and the actual one is obtained.

More generally, if confined plasticity develops at the crack tip, the above decomposition does not hold down to the crack tip. However, excluding the plastic region (or a larger one) around the crack tip allows one to resort to the above decomposition for the elastic displacement. Hence, the above described procedure leads to a unique and consistent definition of the "equivalent elastic crack tip," i.e. the one that would provide the most accurate correspondence between the displacement field of the medium under study and that of a cracked homogeneous elastic medium with the same stress intensity factor.

This strategy straightforwardly applies to any order function. The proper interpretation of this quantity is based on the analysis of a domain where all crack tip non-linearities have been cut out. The elastic description, which contains any negative order (supersingular) field contribution, can thus be extrapolated down to the (elastic equivalent) crack tip. By construction, or definition of the crack tip, the n = -1 or dipolar crack amplitude vanishes.

The quadrupolar term (n = -3) is the first non vanishing intrinsic characteristic of the crack strain field. For all types of non-linearity (e.g. plasticity, damage, cohesive crack tip) the amplitude a -3 exists, and provides a unique and well defined scalar signature of the non-linearities. Higher order terms are associated with fields that vanish more quickly away from the crack tip.

To summarize, for a complete identification of the mechanical state of a material around the crack tip, first a correct position is searched for. Second, stress intensity factors are estimated using the estimated crack tip position.

This procedure may be performed iteratively.

Extraction problem

Let us consider that the displacement field, u, is measured on the surface of a cracked sample. This experimentally measured displacement field is decomposed onto an interpolation basis

u exp (x 1 , x 2 ) = i∈M u exp i n i (x 1 , x 2 ) (13)
In the domain where non-linearities can be neglected, the actual displacement field, u, theoretically belongs to the space generated by the basis functions {φ} as a solution to a Linear Elastic Fracture Mechanics problem. Hence, amplitudes a j i are introduced such that [see Equation [START_REF] Bui | Associated path independent J -integral for separating mixed modes[END_REF]]

u(x 1 , x 2 ) = j a j I φ j I (x 1 , x 2 ) + a j II φ j II (x 1 , x 2 ) (14)
The (best) determination of amplitudes a j i -collectively designated by the vector form {A} -is the objective of the present paper.

The measured displacement field is restricted to the interpolation space generated by the function basis n i , and, in the case of digital image correlation, it is itself obtained by a least squares minimization of an L 2 norm as recalled in the next Section. Introducing the "mass matrix"

M ij = Ω n i (x 1 , x 2 ) • n j (x 1 , x 2 ) dΩ (15)
the measured displacement field reads

u exp (x 1 , x 2 ) = M -1 ij Ω n j (x 1 , x 2 ) • u(x 1 , x 2 ) dΩ n i (x 1 , x 2 ) ( 16
)
where Einstein's convention is used. u exp is thus expressed as

u exp (x 1 , x 2 ) = a j I ϕ j Ii + a j II ϕ j IIi n i (x 1 , x 2 ) (17)
where ϕ j ik can be rewritten in matrix form [Φ] that collects in its rows the components of all φ j i functions projected on the interpolation basis n k . To estimate amplitudes {A}, it is proposed to introduce a set of extracting functions {ψ} in the interpolation space so that

Ω ψ i (x 1 , x 2 ) • Φ jk n k (x 1 , x 2 ) dΩ = δ ij (18)
Numerically, it is more convenient to work with amplitudes Ψ iq of ψ i functions on n q interpolation functions, or

ψ i (x 1 , x 2 ) = Ψ iq n q (x 1 , x 2 ).
The above conjugation property is rephrased as

Ψ iq Ω n q (x 1 , x 2 ) • n k (x 1 , x 2 ) dΩ Φ jk = δ ij (19)
For the set of functions {ψ}, their scalar product with the measured displacement field yields the following identity

ψ i (x 1 , x 2 ) • u exp (x 1 , x 2 ) dΩ = a i (20)
or

{A} = [Ψ][M]{U} (21)
where {U} collects all u exp k . For numerical implementation, it is convenient to compute matrix [V] = [M][Ψ] T , so that the extraction of {A} is given by the following product

{A} = [V] T {U} (22)
For this reason, [V] are referred to as "extractors" or extracting functions.

The present procedure is general and may be tailored to measure as many amplitudes a i as there are degrees of freedom in the measured displacement field. However, in practice, a lack of accuracy will hinder the performance of the extraction. Thus only a limited set of functions φ and ψ will be considered. As a consequence, the conjugation property [START_REF] Rajaram | Application of domain integral methods using tetrahedral elements to the determination of stress intensity factors[END_REF] is not sufficient to compute the extractors. It is proposed to supplement these conditions with another type of constraint, i.e. noise-robustness, as detailed in the following section. A similar approach was also proposed in the framework of the virtual fields method [START_REF] Avril | Sensitivity of the virtual fields method to noisy data[END_REF].

Noise robustness

Being obtained from experimental data (e.g. images), u exp is inevitably corrupted by measurement uncertainties, whose influence on the extracted coefficients is to be reduced as much as possible. A perturbation vector of the component of u exp in the interpolation basis is denoted by {dU}. This perturbation is defined as a zero average noise

{dU} = {0} (23)
with the following correlation kernel

{dU}{dU} T = [C]σ 2 (24)
when the inputs (i.e. the images) of the correlation algorithm are affected by a spatially uncorrelated noise with a standard deviation σ. Section 3 is devoted to the analysis of correlation kernel [C].

The displacement uncertainty subsequently affects the extracted coefficients a i , and by using Equation ( 22), the corresponding perturbation reads

{dA} = [V] T {dU} (25)
Consequently, {dA} is also of zero mean and its variance becomes

{dA}{dA} T = [V] T [C][V]σ 2 (26)
This result quantifies the effect of measurement noise on the extracted coefficients. Noise robustness is achieved through the minimization of the above quantity. For instance, one single amplitude (e.g. stress intensity factor K I ∝ a 1 I ) may be searched for. In that case, only dK I is to be minimized. If several quantities are to be estimated, the minimization of each variance leads to the determination of the corresponding extractor.

Optimal extraction

The procedure of optimal extraction is based on the joint use of the conjugation property [START_REF] Rajaram | Application of domain integral methods using tetrahedral elements to the determination of stress intensity factors[END_REF] and the noise robustness condition [START_REF] Rice | Mechanics of crack tip deformation and extension by fatigue[END_REF]. Let us consider that a 1 is searched for. Only a reduced set of meaningful elastic fields will be considered, namely i ∈ [n 0 ; n 1 ], so that the conjugation conditions are written as

Ψ 1k M kj Φ ij = δ 1i (27)
Reverting to {V 1 } vector, the above equations are recast as

V 1j Φ ij = δ 1i (28)
Under these constraints, the minimum of

R = 1 2 V 1i C ij V j1 (29)
is sought. This constrained minimization problem is solved by using Lagrange multipliers denoted by {Λ}

   [C] [Φ] T [Φ] [0]         {V 1 } {Λ}      =      {0} {∆}      (30)
It is worth noting that if [C] is proportional to the identity matrix, i.e. if the noise is considered as uncorrelated and uniform (white noise), then the optimal extractor is exactly the least squares extractor. Thus least squares extraction is the best extractor for a white noise.

As mentioned in the introduction, the interaction integral [START_REF] Réthoré | Estimation of mixed-mode stress intensity factors using digital image correlation and an interaction integral[END_REF][START_REF] Réthoré | Noise-robust Stress Intensity Factor Determination from Kinematic Field Measurements[END_REF] can be read as a scalar product between the actual displacement and an auxiliary displacement with respect to a given norm involving an interaction Eshelby tensor integrated over the domain affected by a virtual extension of the crack.

In the displacement functional basis, the interaction integral is rewritten in a matrix-vector format as

I int = {U a } T [P]{U} (31)
Thus, if {V int } = [P]{U a }, I int is computed in a similar formalism as {A} and interaction integrals can be recast in the framework proposed herein.

They are associated with a given extracting field v int . Due to the Eshelbian norm in the interaction integral, its associated extracting field is in the subspace of fields that satisfy the extraction conditions. The specific form of the weight function introduces restrictions on the associated extractor field, so that the sub-space to which it belongs is much smaller than the general case considered herein. The optimal version proposed in Ref. [START_REF] Réthoré | Noise-robust Stress Intensity Factor Determination from Kinematic Field Measurements[END_REF] is thus a local minimizer of the noise sensitivity. This property will be illustrated in the example section.

Up to now, the nature of domain Ω was not specified precisely. The only (implicit) property is that it should only contain a region of space where linear elasticity holds. However, it may contain holes, it may not include the crack mouth, it may have a complicated boundary, and the above formulation will still hold. Such a tolerant formulation is suited to the exploitation of image correlation data where information may be lacking in some regions. This is to be contrasted with the interaction integral formalism that requires at least a good capturing of the displacements along the crack faces.

The present framework will be applied to two examples dealing with displacement fields measured by DIC, which is briefly summarized in the next section.

Digital image correlation

In the following, 2D displacement measurements are considered by using pictures f and g of a sample surface. The passive advection of the texture between the reference f and deformed g images reads

f (m) = g(m + u exp (m)) (32)
over the entire pixel map m. From the knowledge of f and g, the measurement problem consists in identifying u exp as accurately as possible. As such, the problem is ill-posed, unless additional assumptions are made on the displacement field so that the information is sufficient to determine u exp with a reasonable accuracy. The quadratic difference is integrated with respect to m over a domain Ω (i.e. a Galerkin approach to the measurement problem)

R cor = Ω [f (m) -g(m + u exp (m))] 2 dΩ ( 33
)
and considered as a residual R cor to be minimized with respect to the set of kinematic unknowns. Adopting a general discretization scheme of the displacement field (see Equation ( 13)) allows for the measurement of the degrees of freedom u exp i .

At this level of generality, many choices can be made to measure 2D displacement fields. When 2D images are considered, 4-noded elements defined on a square grid (i.e. bi-linear functions of x 1 and x 2 [START_REF] Zienkievicz | The Finite Element Method[END_REF]) may be chosen (it is referred to as Q4-DIC [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Châtelier bands[END_REF]). Enriched kinematics (à la X-FEM [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]) may also be considered to describe weak or strong discontinuities associated with Q4 elements (i.e. XQ4-DIC [START_REF] Réthoré | Shear-band capturing using a multiscale extended digital image correlation technique[END_REF][START_REF] Réthoré | Experimental investigation of localized phenomena using Digital Image Correlation[END_REF]). In all the previous cases, no hypotheses are made concerning the underlying material behavior. When the latter is known (e.g. linear and isotropic elasticity), integrated approaches can be developed to directly evaluate mechanically significant quantities (e.g. elastic properties [START_REF] Hild | Digital image correlation: From measurement to identification of elastic properties -A review[END_REF] or even stress intensity factors [START_REF] Roux | Stress intensity factor measurements from digital image correlation: Post-processing and integrated approaches[END_REF][START_REF] Hamam | Stress intensity factor gauging by digital image correlation: Application in cyclic fatigue[END_REF][START_REF] Réthoré | An extended and integrated digital image correlation technique applied to the analysis fractured samples[END_REF]). This last route is not followed hereafter, rather so-called Q4 and XQ4 kinematics are measured and subsequently post-processed.

The interpolation functions being chosen, the minimization of R cor is carried out using a modified Newton procedure

[N]{U i } = {b i } (34)
that is iteratively solved until convergence is reached. The following definition holds for the search direction [N] of the modified Newton procedure and the residual {b i } at iteration i

N kl = Ω (n k (m) • ∇f (m)) (n l (m) • ∇f (m)) dΩ (35)
and

b i k = Ω (n k (m) • ∇f (m)) f (m) -g(m + (u exp (m)) i ) dΩ ( 36
)
The gradient of the image as well as sub-pixel interpolations have been the subject of publications in the literature [START_REF] Sutton | Effects of subpixel image restoration on digital correlation error estimates[END_REF][START_REF] Schreier | Systematic errors in digital image correlation caused by intensity interpolation[END_REF]. A standard choice is a finite difference gradient and cubic spline interpolation.

Let us now consider that the images are affected by noise. It is assumed that f is a pure (noise free) reference image and noise is reported on fg. This noise is considered as spatially uncorrelated, of zero mean and its variance is equal to twice the variance of the sensor noise σ 2 c . Typical values for σ c are 1 or 2 gray levels for 8-bit pictures. As a consequence, the image noise is reported on the LHS {b i } of Equation (34) for which the variance

reads b i k b i l = 2σ 2 c Ω (n k (m) • ∇f (m)) (n l (m) • ∇f (m)) dΩ (37)
The search direction [N] then leads to the displacement perturbation

{dU}{dU} T = 2σ 2 c [N] -1 [N][N] -1 (38)
The correlation kernel of the displacement uncertainty is identified as

[C]σ 2 = 2σ 2 c [N] -1 (39)
Covariance matrix [C] is directly related to the image texture (i.e. its gradient) and the chosen measurement basis n k . The joint use of the last results (Equations ( 39) and ( 30)) allows for the computation of the optimal extractor precisely suited to a specific image, and a chosen kinematic basis.

Examples

The present section shows the application of the above formalism to two examples. In the first one -a fatigue test on a steel specimen -the above described procedure is used to identify the crack tip position and stress intensity factors. The second one deals with a brittle ceramic material for which a very low signal-to-noise ratio is observed. In that case, the extraction of the mode I stress intensity factor is addressed, and a comparison is performed with alternative techniques.

Fatigue test on a steel CCT specimen

This first example is devoted to the analysis of a CCT specimen made of XC48 (or C45) steel with Young's modulus of 190 GPa and Poisson's ratio of 0.3. The cyclic yield strength σ y (offset: 0.2 %) is equal to 210 MPa.

The specimen is subjected to cyclic tension with a load ratio R = 0.4. In the present analysis, only the stage corresponding to the maximum load level is considered after about 300,000 cycles for which the crack size 2a = 14.5 mm. The sample thickness is equal to 5 mm. The vertical and horizontal components of the displacement field obtained by using 12 × 12-pixel Q4 elements are shown in Figure 1. The conversion factor p is 2.08 µm/pixel.

In this example, a process zone is expected. The exact position of the crack tip along the crack mouth is searched for. This location is critical if a stress intensity factor is to be evaluated [START_REF] Roux | Stress intensity factor measurements from digital image correlation: Post-processing and integrated approaches[END_REF]. Thus the supersingular component K -1

I

has to be extracted. It allows for an estimate of the crack tip location [START_REF] Hamam | Stress intensity factor gauging by digital image correlation: Application in cyclic fatigue[END_REF][START_REF] Réthoré | An extended and integrated digital image correlation technique applied to the analysis fractured samples[END_REF][START_REF] Henninger | Enriched kinematic fields of cracked structures[END_REF]. Further, a mask is adopted to account for linear elastic displacement fields, i.e. a 20-pixel wide band and a core 20 pixels in radius are omitted for the extraction. In the sequel, the domain size is the diameter of the disk that defines the outer boundary of the mask.

Crack tip location

In Figure 2(a), the estimation of the shift x I with respect to the prescribed position of the crack tip along the crack faces is plotted. A linear dependence is obtained with a -1 slope. From Equation ( 12), the ratio between the shift of the crack tip x I and the crack tip position used to compute K & M potentials is theoretically equal to -1. This is confirmed by Figure 2(a) for both optimal and least squares extraction techniques, the latter being 1 pixel or about 2 µm higher. The variation of the crack tip shift (and thus of K -1 I ) is evaluated by using the results of Figure 2 As observed for the stress intensity factor, the optimal extraction allows for a reduction of this noise sensitivity compared with the results obtained with the least squares approach. Considering that the noise sensitivity of K I has only a small contribution to the noise sensitivity of the identification of the crack tip position, the uncertainty level of the shift estimation is about 0.03 µm/gray level for the optimal extraction.

SIF extraction

The results on the stress intensity factor extraction and their noise sensitivity are presented in Figure 3. In the present case, the signal-to-noise ratio is high. Consequently, least squares and optimal extractions give very similar results for K I and K II , namely, a value of 20.7 MPa √ m is obtained for K I , and 0.05 MPa √ m for K II . Because of the high signal-to-noise ratio, a very low dependence to the size of the integration domain is observed as well as very low noise sensitivity. A 0.0011 MPa √ m/gray level standard uncertainty is obtained for the least squares extraction, and 0.0009 MPa √ m/gray level for the optimal approach. The region of interest covers an area of 706 × 646 pixels approximately centered about the crack tip. The displacement field is measured by using 16 × 16-pixel Q4 elements with a discontinuous enrichment of the kinematic basis [START_REF] Réthoré | From pictures to extended finite elements: Extended digital image correlation (X-DIC)[END_REF]. Figure 6 shows the two components of the measured displacement field. Due to the poor quality of the image (the rough surface of the material is used as a random texture for image correlation), and also the small magnitude of the displacement fluctuations, a low signal-to-noise ratio is observed. A normal displacement discontinuity of about 0.5 pixel (at most) is captured at the bottom of the region of interest (i.e. the crack is invisible).

Extracting fields

The conversion factor p is 1.85 µm/pixel.

Extraction analysis and 3D effects

In this study, only the optimal extractor with an L 2 projection is used.

A mask with a disk of radius r 0 = 20 pixels, in addition to a 20 pixel wide masked zone around the crack surface is adopted. Crack tip shift and process zone width are shown in Figure 7. Compared to the results of the previous section, the extracted quantities suffer from fluctuations of higher magnitude when the domain size is increasing. A value of the crack tip shift with respect to its initial position of about -67 µm is obtained. When the change of the crack tip shift with respect to the assumed crack tip position is observed, the expected linear dependence is rather well retrieved.

Different arguments can be invoked in order to explain these unexpected results. The first point is that the displacement field is measured on a free surface, and hence the analytic 2D plane stress displacement field that is used is not the actual one. A second point is that a post-mortem observation of the specimen and the crack surface (Figure 8) shows that the crack front is not perfectly orthogonal to the free surfaces on which the displacement was measured in the correlation analysis. Thus, an additional 3D effect is expected. To evaluate its influence, a 3D elastic analysis is carried out on the entire specimen with idealized crack geometry and boundary conditions. A planar crack is considered and the front is modeled as a parabola such that the crack length on the free surface is half of the specimen width and 180 µm longer in the middle section, consistent with post-mortem observations (Figure 8). A straight front (normal to the observation face) is also considered as a reference case.

The results are shown in Figures 9 and10. In Figure 9, the displacement field on the free surface for a straight crack front is depicted. For this analysis, the element size in the region of interest is equal to 16 × 16 pixels (for a pixel size of 2.11 µm). In Figure 10(b), the same quantities as for Figure 7(b) are given as functions of the domain size. Three cases are considered. First, a straight crack front is used and the extraction analysis is carried out on the free surface and also on the mid-plane of the specimen (under plane strain conditions). Second, the displacement field on the free surface is analyzed for the parabolic crack front that fits the actual geometry observed on the postmortem photograph. In Figure 10(a), stress intensity factors are normalized by the value for the mid-plane and the straight front. It is observed that the crack front geometry has a strong influence not only on the mean value of K I but also on its dependence with the domain size.

Further, opposite tendencies are obtained for the three different cases, namely, for the analysis on the mid-plane of a crack with a straight front, K I is almost constant and is chosen as a reference value K Io , on the free surface the value of K I increases from 0.95 to 0.97 the reference value K Io whereas for the curved front K I decreases from 1.11 to 1.09. The same trend is obtained for the shift and the process zone width. A 0 shift is obtained for the straight front for the examined section but a shift of about -20 µm is obtained for the parabolic crack front (Figure 10(c)), meaning that the depth of influence of the front geometry is about 350 µm (or 165 pixels), which is about 0.12 times the specimen thickness, and 0.07 times its width for an apparent crack length of 0.5 times its width.

Even though the previous analysis may only be qualitative, it is proven that 3D effects may exhibit supersingular contributions. Compared with the results obtained in the previous paragraph for the experimental specimen, the same trends are obtained, namely, the crack tip position is estimated in front of its emergence on the free surface due to curvature effects and a process zone size of about 40 to 50 µm is estimated.

A mask is used to select the domain over which displacements are used in the analysis. The outer boundary of the mask is now chosen to be the largest disk included inside the region of interest. In addition, the immediate vicinity of the crack faces is excluded from the analysis as well as a small disk of radius r 0 centered on a rough prior determination of the crack tip.

Figure 11 shows the influence of this inner radius r 0 . It is observed that this excluded crack tip region has a significant influence as could have been anticipated by the singular nature of the fields at the crack tip. However, remarkably, for r 0 less than about 175 µm, no systematic trend is observed (although fluctuations are present).

SIF extraction

In this last paragraph, the gain brought by an estimation not only of the displacement field but also of its discontinuity is evaluated. The analysis is carried out using the crack tip position determined in the previous The RMS error for an image noise of 2 gray levels is displayed with error bars in Figure 12. For the largest integration domain, the uncertainty reduces to 0.02 MPa √ m for the optimal extraction whereas it is equal to 0.05 MPa √ m when the least squares technique is applied.

These results are compared to those obtained in Ref. [START_REF] Réthoré | Noise-robust Stress Intensity Factor Determination from Kinematic Field Measurements[END_REF] when using a standard or optimized interaction integral. Figure 13(a) shows the change of K I for four extraction techniques, namely, least squares and optimal techniques, interaction integral with a conical virtual crack extension field (VCE)

and interaction integral with an optimal VCE. Although the interaction integral is theoretically domain independent [START_REF] Parks | A stiffness derivative finite element technique for determination of crack tip stress intensity factors[END_REF], larger variations with the domain size are observed due to a higher sensitivity to noise. Yet, the values obtained with the four techniques are consistent for large integration domains, namely, 2.9 ± 0.09 MPa √ m and 3.0 ± 0.06 MPa √ m for the interaction integral with conical and optimal VCEs. of the integration domain) but they differ from the optimal interaction integral D -1.5 and the standard interaction integral D -1 . In Ref. [START_REF] Réthoré | Noise-robust Stress Intensity Factor Determination from Kinematic Field Measurements[END_REF], the D -1.5

decrease is shown to be the optimal rate for the interaction integral. Let us note that the formalism proposed herein allows for higher reduction rates of noise sensitivity (i.e. D -2 ). This is due to the fact that differential operators are involved in the interaction integral and reduce the noise reduction rate of any derived extractor. The difference between least squares and optimal extraction appears to be a multiplicative constant only. Due to the very brittle nature of the ceramic, one may think of relaxing the orthogonality conditions with respect to supersingular modes. As a consequence, the optimal field becomes more uniform allowing all displacement degrees of freedom to contribute to the determination of the displacement field. They are used to obtain the crack tip position.

Two experimental applications have been proposed. In the first one, the small thickness of the specimen legitimates the use of a two-dimensional description, and the results are shown to improve the robustness and accuracy of the measurements of fracture properties. In the second example, although the large scale results seem consistent, a more pronounced domain dependence still remains, which is attributed to three dimensional effects. In fact, this difficulty is a major one, as far as the measurement itself relies on a description that is not exact. The 2D analysis is however used here in order to conform with the standard practice, yet a better suited field (as well as an amplitude equivalent to a SIF that can be quantitatively used) should be more thoroughly investigated. Examining domain integral approaches for stress intensity factor extraction within the framework proposed herein, unexpected results are obtained.

Re-read as an extracting field, the interaction integral extracts stress intensity factors predominantly using displacement data at the nodes located close to the crack faces. In the present examples, it is emphasized by the use of a conical virtual crack extension fields. More regular virtual crack extension fields could be used [START_REF] Rajaram | Application of domain integral methods using tetrahedral elements to the determination of stress intensity factors[END_REF] to enlarge the useful data area around the crack tip, yet concentrated weights still remain along the crack faces.

Last, it is to be emphasized that the general procedure presented herein to compute optimal extractors is also suited to dealing not only with uncertain data, but also with kinematic fields being defined on arbitrary shape and topology domains, a property that is not found for, say, interaction integrals.

- 1 I 3 I{U} 1 , u 2 1 } 1 z

 131211 Amplitude of a displacement field decomposed over Williams' series {A} Vector representation of a n i A T Transpose of matrix or vector A {dA} Vector representation of perturbation on a n i {b i } Vector form of the correlation residual at iteration i Derivation with respect to k th spatial coordinate K i Stress intensity factor in mode i = I or II K Io Reference value for K I δK Perturbation on the stress intensity factor K K First mode I super-singular factor K -Third mode I super-singular factor Vector representation of the nodal displacement field {dU} Vector representation of the nodal displacement perturbation u Displacement field u a Auxiliary displacement field {U a } Vector representation of the nodal auxiliary displacement field u exp , u exp i Experimental displacement field, and i-th component u Displacement jumps in the tangential and normal directions, respectively [V] Matrix representation of the nodal extractor field {V Vector representation of the nodal extractor field associated with a 1 {v int } Extractor vector field associated with I int {V int } Vector representation of the nodal extractor field associated with I intx iLocal coordinate system at the crack tip (i = 1, 2)x I Shift of the tip from its actual position along x Position on the interpolation basis n k[Φ] matrix that collects in its rows the components of all φ j i functions projected on the interpolation basis n k ψ i Extracting functions {Ψ i } Vector representation of the nodal extracting function θ Polar angle of the current point with respect to the crack tip axis τ K Noise sensitivity of mode I stress intensity factor τ K -1 Noise sensitivity of mode I first super singular factor

Section 4 ,2. Problem formulation 2 . 1 .

 421 two examples are discussed. The first one is dedicated to a fatigue crack in a steel sample where a confined plastic zone exists at the crack tip. Stress intensity factors and crack tip position are estimated using the optimal extraction of appropriate crack displacement field parameters. Second, the extraction of stress intensity factors is examined from a low signal-to-noise ratio displacement field of a brittle material. Crack displacement fields Let us consider a semi-infinite planar crack whose normal is aligned along the x 2 axis, and crack front parallel to the x 3 axis. The crack corresponds to the domain x 2 = 0 and x 1 < 0. Crack displacement fields u in the linear elastic regime are computed as a double infinite series of Kolossov-Muskhelishvili (K & M) functions [31] φ n I and φ n II for modes I and II. The complex plane formulation reveals convenient in the present case. A point z in the complex plane x 1 + ix 2 is represented in polar coordinates (r, θ)

  All these fields are self-similar, and thus are homogeneous functions of the distance to the crack tip r. It is convenient to label them with the integer index n that indicates the exponent of the power-law dependence of the displacement with r, that is n/2. With this convention, all elastic displacement fields are generated with -∞ < n < ∞. All odd indices label fields with a discontinuity across the crack mouth, whereas even indices correspond to continuous fields. To mention just a few examples, n = 0 corresponds to translations either parallel (hence termed I) or perpendicular (labeled II) to the crack path, n = 1 are the classical mode I and mode II displacement fields whose amplitudes are proportional to the corresponding stress intensity factors. The case n = 2 gives access to either the T-stress component for mode I, or rigid body rotation about the crack tip for mode II

  (b). A typical variation of 5 µm is obtained for both extraction techniques. The range within which K -1 I is less dependent on the domain size appears lower than the independence range of the stress intensity factors (see Figure 3(a)), which could have been expected considering the higher variation of the first supersingular field compared with the usual √ r-dependence of the near crack tip field.

Figure 2 (

 2 Figure 2(c) shows the noise sensitivity of the extraction of K -1 I . As observed for the stress intensity factor, the optimal extraction allows for a

Figures 4 - 1 I 4 . 2 .

 4142  show the components of the extracting fields for K 1 I and K -1 I using least squares and optimal extractions. As for K I , the optimal fields spread the weight allowed to each displacement degree of freedom depending on its associated level of confidence. First, due to the divergence of supersingular fields, n = -1, in the vicinity of the crack tip, their associated extracting fields weight more intensively the displacement data around the core of the mask. Second, the angular dependence, bipolar for n = -1, could have been derived from the expression of the supersingular fields. Compared

  paragraph. First, orthogonality conditions with K & M "harmonics" in the interval [-3; 5] are prescribed. The mode I stress intensity factors measured using least squares technique and the proposed optimal extraction are shown in Figure 12. K I is plotted as a function of the diameter of the circular integration domain around the crack tip. A value of 3.1 MPa √ m is obtained.

Figure 13 (

 13 Figure 13(b) shows the standard uncertainty for K I with respect to the domain size. Power law decreases are obtained for the four techniques. The least squares and the optimal extractions vary as D -2 (D being the diameter

Figures 14 and 15

 15 Figures 14 and 15 show the extracting fields of mode I stress intensity factors for the four above mentioned techniques. The least squares extracting field (Figure 14(b)) gives a considerable weight to the displacement values very close to the crack tip.The optimal version spreads the weight over a larger neighborhood, and accounts also for the non-homogeneous noise sensitivity of the displacement field measurement (Figure14(a)). This effect allows for a lower noise sensitivity of the extracted stress intensity factor K I . By looking at the extracting fields related to the interaction integrals, a considerable weight is carried by the displacement values not only close to the crack tip but also along the crack faces (Figure15(b)). Further, a flower shape around the crack tip is observed. Because of the presence of differential operators in the integration integral, the shape of the extracting field is more

Figure 16 (

 16 Figure 16(a) shows the corresponding extracting field. Concerning the least squares extracting field, Figure 16(b) also shows less concentrated weighting but its noise sensitivity increases significantly (Figure 17(b)). Consequently, K I values are not consistent with previous determinations (Figure 17(a)), whereas the optimal extraction with reduced orthogonality constrains gives consistent K I values of 3.1 MPa √ m with a very low uncertainty level, namely 0.01 MPa √ m. This is presumably another illustration of the above discussed 3D effects.

  The performances of the proposed optimal extraction are illustrated and compared to existing techniques through two examples. A significant improvement of the noise sensitivity is obtained by reaching a global minimizer of the effect of noise instead of a local one (for least squares or interaction integral techniques). Domain independence is obtained for stress intensity factors in spite of the mechanically non-admissible nature of measured displacement fields i.e. they do not satisfy the balance of momentum equation. Domain independence is also obtained for supersingular components but over a reduced domain size range.

Figure 1 :Figure 2 :

 12 Figure 1: Vertical (a) and horizontal (b) components of measured displacement field by Q4-DIC for the CCT specimen, expressed in pixels (1 pixel ↔ 2.08 µm).

Figure 3 :

 3 Figure 3: Stress intensity factors (a) and their noise sensitivity (b) for the CCT specimen obtained with a least squares and optimal extraction performed with an L 2 projection and orthogonality conditions with K & M harmonics in [-3; 5].

Figure 4 :Figure 5 : 1 IFigure 7 :

 4517 Figure 4: Vertical (top) and horizontal (bottom) components of the extracting field normalized by 2µ √ pπ obtained for an optimal extraction (left) a least squares extraction (right) of K I performed with an L 2 projection and orthogonality conditions with K & M harmonics in [-3; 5].

Figure 8 :

 8 Figure 8: (left) Reference picture of the analyzed region of interest in the SiC sample (the fine notch can be seen on the bottom side). (right) post-mortem view of the crack surface.The pre-notch, the fine notch, and the precrack surfaces are clearly observable. Note the curvature of the precrack front.

Figure 9 :

 9 Figure 9: Vertical (a) and horizontal (b) components of the displacement field on the free surface of a three-dimensional elastic solid expressed in pixels obtained with a straight crack (1 pixel ↔ 2.11 µm).

Figure 10 :

 10 Figure 10: Normalized stress intensity factors (a) and shift (b) of the 3D elastic solid for a parabolic crack front, a straight front on the free surface and in the center section of the specimen (1 pixel ↔ 2.11 µm). The crack tip shift is also plotted as a fonction of the crack tip position (c) (a dotted line of slope -1 is shown as a guide to the eye).

Figure 13 :Figure 14 :Figure 15 :

 131415 Figure 13: Comparison of mode I stress intensity factors and noise sensitivities for the SiC specimen obtained with different extraction techniques including previous results with the interaction integral [24].

Acknowledgments

This work was funded by a CETIM Foundation grant entitled PROPA-VANFIS: Advanced methods for the experimental and numerical analyses of crack propagations under complex loadings, and by Agence Nationale de la Recherche through grant RUPXCUBE (ANR-09-BLAN-0009-01).


Baldrik Faure
email: baldrik.faure@inrets.fr

Olivier Chiello
email: olivier.chiello@inrets.fr

Marie-Agnès Pallas
email: marie-agnes.pallas@inrets.fr

Christine Serviere
email: christine.serviere@gipsa-lab.grenoble-inp.fr

Christine Servière

Caractérisation du rayonnement acoustique d'un rail à l'aide d'un réseau de microphones -Simulations

come

Introduction

Dans un contexte de développement des transports ferroviaires urbains et péri-urbains, la réduction des nuisances sonores apparaît comme un enjeu majeur. Afin d'agir efficacement à la source, il est indispensable d'identifier et d'étudier les sources responsables de ces nuisances. Pour des vitesses inférieures à 300 km/h, le bruit de roulement constitue la source principale du bruit ferroviaire [START_REF] Mellet | High speed train noise emission : Latest investigation of the aerodynamic/rolling noise contribution[END_REF].

La modélisation du bruit de roulement a été initiée par Remington et reprise par Thompson [START_REF] Thompson | A review of the modelling of wheel/rail noise generation[END_REF], notamment à travers le modèle TWINS [START_REF] Thompson | Experimental validation of the TWINS prediction program for rolling noise, Part 1 : description of the model and method[END_REF]. Dans ce modèle, le déplacement vertical relatif au niveau du contact roue/rail (dû aux irrégularités des surfaces de roulement) engendre des vibrations qui se propagent dans la voie. Le rayonnement acoustique des éléments tels que le rail, les roues et les traverses constitue le bruit de roulement. Ce modèle a montré que la contribution du rail au bruit de roulement reste importante pour des fréquences comprises entre 500 et 1200 Hz [START_REF] Thompson | A review of the modelling of wheel/rail noise generation[END_REF]. Ces prédictions ont en partie pu être observées lors de campagnes de mesures utilisant des méthodes dérivées du principe de formation de voies ; une source importante et étendue a été identifiée au niveau du rail [START_REF] Pallas | Tram noise emission : spectral analysis of the noise source contributions[END_REF]. Cependant, le principe du traitement d'antenne tel qu'il est utilisé actuellement est mal adapté à ce type de source [START_REF] Kitagawa | Comparison of wheel/rail noise radiation on Japanese railways using the TWINS model and microphone array measurements[END_REF]. En effet, le comportement vibro-acoustique du rail, en particulier la propagation de vibrations sur une très grande distance dans certaines gammes de fréquences, conduit à le considérer comme une source étendue ayant une grande longueur de cohérence.

Dans l'étude présentée ici, le rayonnement acous-tique du rail est caractérisé par des paramètres vibratoires (amplitudes des excitations, nombres d'ondes complexes) que l'on cherche à estimer à l'aide d'un réseau de microphones. La méthode proposée consiste à minimiser l'erreur entre les matrices interspectrales modélisée et mesurée sur les microphones de l'antenne. Dans un premier temps, le modèle vibro-acoustique retenu pour le rail est présenté. Les caractéristiques du champ acoustique qui en résultent sont ensuite étudiées. Enfin, par minimisation de l'erreur entre modèle et mesures, certains paramètres vibratoires sont évalués dans le cas simple d'une ou plusieurs excitations fixes, en ne considérant qu'un seul type d'onde vibratoire. Des résultats de simulations viennent appuyer et compléter l'étude analytique.

Modélisation du rail 2.1 Modèle vibratoire

Un modèle continu classique est adopté [START_REF] Knothe | Modeling of railway track and vehicle/track interaction at high frequencies[END_REF] dans lequel le rail est assimilé à une poutre épaisse de Timoshenko. Pour cette étude, une pose sur ballast est retenue. Les systèmes ballast-traverse-semelle sont modélisés par des impédances mécaniques localisées de type ressort-masse-ressort notées s(ω). Pour des longueurs d'onde inférieures à l'espacement a entre traverses (jusqu'à environ 1500 Hz), ces supports discrets peuvent être remplacés par un support continu de raideur linéique équivalente s ′ (ω) = s(ω)/a. Une représentation schématique de ce modèle est donnée figure 1. Le déplacement vertical du rail en réponse à une force unitaire en z 0 , appelé fonction de Green, s'écrit [START_REF] Hamet | Railway noise : use of the Timoshenko model in rail vibration studies[END_REF] :

G ω (z, z 0 ) = F d e -k d |z-z0| + iF p e -ikp|z-z0| (1)
Dans l'équation La figure 2 représente l'atténuation en dB/m de ces deux ondes ainsi que les fréquences caractéristiques du système. Ces fréquences correspondent aux résonances des différents éléments de la voie. Pour une onde de type e -γ|z-z0| , l'atténuation ∆ en dB/m est calculée ainsi : ∆ = ℜ(γ) × 20 log e. Avec ℜ(γ) la partie réelle de γ.

Modélisation acoustique

Comme dans la référence [START_REF] Thompson | Investigation into the validity of two-dimensional models for sound radiation from waves in rails[END_REF], le rail est assimilé à une ligne de monopôles répartis continûment sur son axe comme le montre la figure 3. La pression acoustique en un point de l'espace de coordonnées cylindriques (r, z) est alors donnée par l'intégrale :

p(r, z) = ∫ +∞ zs=-∞ iωdQ(z s ) e -ikrs(z,zs) 4πr s (z, z s ) (2)
avec : z s l'abscisse d'un monopôle élémentaire sur l'axe z, r s la distance entre le monopôle d'abscisse z s et le point d'observation. dQ(z s ) est le débit massique élémentaire de la source placée à l'abscisse z s ; il est proportionnel à la vitesse vibratoire du rail en ce point. Pour effectuer le calcul numérique de la pression, l'intégrale (2) est discrétisée. La distance D zs entre deux monopôles successifs doit être suffisamment petite par rapport aux longueurs d'onde vibratoires dans le rail et dans l'air. Une précision suffisante est assurée pour :

D zs = 1 5 min(λ air , λ rail ) (3)
en imposant D zs < 0.1 m pour assurer une bonne estimation proche du rail.

3 Particularités du champ acoustique rayonné par le rail La figure 2 montre que l'atténuation de l'onde propagative varie avec la fréquence de l'excitation. L'onde de champ proche est quant à elle fortement atténuée quelle que soit la fréquence. Les cartographies suivantes illustrent les conséquences de ces différents comportements sur le champ acoustique rayonné par le rail pour deux fréquences caractéristiques.

Cas simple : un seul contact d'amplitude inconnue

On considère dans un premier temps le cas d'un contact unique à l'abscisse z 0 . Les nombres d'onde complexes dans le rail sont connus pour l'onde propagative et l'onde de champ proche ; les paramètres dynamiques de la voie ont par exemple été déterminés au marteau de choc. Seule l'amplitude du champ, liée à l'amplitude de l'excitation et au facteur de rayonnement, est inconnue. Pour ce cas simple, une écriture analytique de la solution du problème de minimisation est possible. Le champ de pression p(r, z) donné par l'équation (2) est proportionnel à un champ de pression élémentaire p 0 (r, z) :

p(r, z) = A • p 0 (r, z) = A • ∫ +∞ zs=-∞ -ω 2 G ω (z 0 , z s ) e -ikrs (z,zs) 4πrs(z,zs) d zs (5)
où A est l'amplitude complexe du champ élémentaire. Le critère Λ peut donc s'écrire :

Λ = ∑ m,n |Γ m,n -|A| 2 (P 0 • P 0 † ) m,n | 2 = ∑ m,n |Γ m,n -|A| 2 M m,n | 2 (6)
avec : P 0 le vecteur des N mic pressions élémentaires sur l'antenne M la matrice spectrale élémentaire modélisée. En posant |A| 2 = α (α > 0), le critère (6) atteint un minimum lorsque sa dérivée par rapport à α est nulle :

∂Λ ∂α = 0 ⇔ α = ∑ m,n M m,n Γ * m,n ∑ m,n |M m,n | 2 (7)
où Γ * m,n est le conjugué de Γ m,n . Remarque : A est supposé ici déterministe. Dans le cas d'une variable aléatoire, les expressions ( 6) et [START_REF] Hamet | Railway noise : use of the Timoshenko model in rail vibration studies[END_REF] restent valables : il suffit de substituer à |A| 2 la variance σ a 2 de A.

Performances de l'estimateur en présence de bruit

Dans le cas où un bruit additif spatialement blanc, centré et de puissance moyenne σ 2 , se superpose aux signaux mesurés, la matrice spectrale est donnée par :

Γ = Γ + σ 2 I Nmic (8)
avec I Nmic la matrice identité de dimension N mic (nombre de microphones).

En pratique, Γ est estimée par une moyenne sur K observations. Soit α l'estimée de α en présence de bruit. En utilisant les équations [START_REF] Hamet | Railway noise : use of the Timoshenko model in rail vibration studies[END_REF] et [START_REF] Thompson | Investigation into the validity of two-dimensional models for sound radiation from waves in rails[END_REF], on peut écrire le biais relatif (équation 9) et la variance relative (équation 10) de l'estimateur :

α -α α = 1 N mic • 1 R s/b (9)
V ar(α)

α 2 = 1 K • [ 2 N mic R s/b + 1 N mic 2 R s/b 2 ] (10)
où R s/b est le rapport signal à bruit linéaire, défini par le rapport de σ 2 et de la puissance moyenne sur les microphones de l'antenne.

Les équations ( 9) et (10) montrent que les performances de l'estimateur ne dépendent pas de la fréquence. Même lorsqu'il y a autant de bruit que de signal, l'erreur d'estimation sur α reste inférieure à 1 dB en moyenne.

Sous cette hypothèse de bruit spatialement blanc, seule la diagonale des matrices spectrales porte la contribution du bruit, introduisant de ce fait le biais évoqué précédemment. Pour pallier à ce problème, il est possible de modifier le critère à minimiser, en excluant les éléments diagonaux des matrices spectrales. L'estimateur ainsi défini n'est plus biaisé, mais une attention toute particulière devra être portée aux résultats obtenus dans les cas plus complexes.

Performances de l'estimateur en présence d'erreurs sur le modèle

La méthode décrite précédemment repose sur la minimisation d'une erreur entre mesures et prévisions issues d'un modèle. Jusqu'à présent, les paramètres de ce modèle sont supposés connus, à l'exception de l'amplitude A de l'excitation. Cette partie présente les résultats d'une étude de robustesse vis-à-vis d'erreurs sur ces paramètres, sans bruit de mesure sur les microphones.

La matrice spectrale mesurée est simulée à partir du modèle vibro-acoustique. Pour la matrice spectrale modélisée, des erreurs sont introduites sur les paramètres suivants : position du contact, atténuation ou nombre d'onde de l'onde propagative. Dans les illustrations numériques, pour une fréquence f donnée l'antenne comporte 13 microphones, espacés de façon à vérifier le critère de Shannon spatial (d ≤ λ air /2) ; sa longueur L est donc spécifique à chaque fréquence. L'antenne est disposée parallèlement au rail à 3 m de celui-ci. La position du centre de l'antenne z c sur cet axe est un paramètre qui varie de 0 à 5 m (l'origine z c = 0 étant face au contact z 0 ).

Lorsqu'une erreur est introduite sur l'atténuation de l'onde propagative, la robustesse est vérifiée pour les fortes atténuations : pour une erreur de ±10 dB/m sur l'atténuation, on constate moins de 1 dB d'erreur sur l'estimation de α. Pour les faibles atténuations, une imprécision sur l'atténuation peut se traduire par une forte erreur sur l'estimation de α, et ceci d'autant plus que l'antenne est excentrée du point de contact (figure 6). Plus généralement, c'est lorsque l'antenne est en face du contact que l'erreur d'estimation est minimale. résultats montrent que la position z 0 du contact est un paramètre à connaître avec précision, surtout pour les basses fréquences. En effet, en hautes fréquences, lorsque les ondes rayonnées par le rail sont planes, la position de l'antenne dans un tel champ invariant dans l'espace importe peu. En pratique, cette information de position du contact par rapport à l'antenne sera a priori assez bien connue.

Amplitude et atténuation inconnues

Étude du critère à minimiser

On suppose à présent que deux paramètres sont inconnus, à savoir l'amplitude de l'excitation et l'atténuation de l'onde propagative dans le rail. La solution de ce problème, par une approche de minimisation de l'erreur entre les matrices spectrales mesurée et modélisée, n'a pas de solution analytique simple. La résolution nécessite l'utilisation d'une méthode d'optimisation. Au préalable, nous analysons le comportement du critère à minimiser, et notamment l'existence éventuelle de minima locaux ou de zones plus vastes où la fonction coût varie faiblement autour de la solution. On considère une antenne microphonique de structure similaire à celle de la section 4.3.

Quelle que soit la fréquence observée, et quelle que soit la position z c du centre de l'antenne par rapport au contact, on constate l'existence d'une vallée le long de laquelle Λ varie très peu. La figure 8 tion que nous utilisons dans un premier temps, est un algorithme de programmation séquentielle quadratique (méthode SQP), particulièrement adapté aux problèmes non linéaires sous contraintes comme le nôtre [START_REF] Biggs | Constrained Minimization Using Recursive Quadratic Programming[END_REF]. Même lorsque l'on initialise à proximité du minimum, l'algorithme converge vers une mauvaise solution quelle que soit le cas étudié (fréquence, position de l'antenne). Le problème ainsi défini est mal conditionné.

N contacts d'amplitude inconnue

On considère à présent N contacts, associés à N roues sur le rail. Les contacts sont supposés décorrélés et d'amplitudes inconnues. Par extension de la section 4.1, chaque contact d'indice u à l'abscisse z u (u ∈ [[1, N ]]) engendre un champ acoustique qui s'écrit sous forme d'un champ élémentaire p u (r, z) multiplié par une amplitude complexe A u . Le champ acoustique total rayonné par le rail excité par l'ensemble de ces N contacts s'écrit :

p(r, z) = ∑ N u=1 A u • p u (r, z) = ∑ N u=1 A u • ∫ +∞ -∞ -ω 2 G ω (z u , z s ) e -ikrs (z,zs) 4πrs(z,zs) d zs (12)
Sous l'hypothèse de décorrélation des N excitations, la matrice spectrale modélisée sur les N mic capteurs s'écrit :

Γ mod = N ∑ u=1 σ u 2 Γ u (13)
avec : σ u 2 la variance de l'amplitude du contact u, Γ u la matrice spectrale sur les N mic capteurs relative au contact u.

Des équations (13) et (4) on déduit le critère Λ pour N contacts :

Λ = ∑ n,m |Γ n,m | 2 -2 ∑ u α u U u + ∑ u,v α u α v V u,v (14) avec
: U u = ∑ m,n (P u ) * m Γ m,n (P u
) n , P u le vecteur de la pression unitaire modélisée sur l'antenne, relatif au contact u,

V u,v = ∑ m (P u ) m (P v ) * m 2 , α u = σ u 2
Ainsi, d'après l'équation (14), trouver les α u qui minimisent Λ revient à résoudre :

∀u ∈ [[1, N ]] ∂Λ ∂α u = 0 (15)
D'après les équations (14) et (15) minimiser le critère Λ pour N contacts revient à résoudre le système matriciel :

   V 1,1 • • • V 1,N . . . . . . . . . V N,1 • • • V N,N       α 1 . . . α N    =    U 1 . . . U N    (16)
La solution d'un tel système est unique si et seulement si la matrice V est inversible. Dans ce cas, si les α u trouvés sont positifs, ils sont également solution du problème d'optimisation. Dans les autres cas, lorsque la matrice V n'est pas inversible ou que les α u trouvés ne sont pas tous positifs, il faudra utiliser un algorithme d'optimisation pour minimiser le critère (14) sous contrainte de positivité des α u . En l'absence de bruit, les résultats obtenus par inversion de la matrice V sont bien les solutions du problème d'optimisation, sur l'ensemble de la gamme de fréquence étudiée (jusqu'à 3000 Hz).

Conclusion

Le caractère étendu du rail en fait une source acoustique difficile à caractériser avec les méthodes de traitement d'antenne classiques ; le champ acoustique qu'il rayonne présente en effet des propriétés trop hétérogènes avec la fréquence. Dans cet article, une méthode d'optimisation paramétrique basée sur la minimisation d'un critère d'erreur quadratique est proposée ; il s'agit de l'erreur commise entre la matrice spectrale mesurée sur l'antenne et une matrice spectrale calculée à partir d'une modélisation vibro-acoustique du rail. Des simulations ont permis de vérifier la robustesse de cette méthode aux erreurs de modélisation pour un cas simple où seule l'amplitude de l'excitation est inconnue. Lorsque deux paramètres sont inconnus, l'utilisation de plusieurs positions pour l'antenne rend le problème d'optimisation mieux conditionné. Enfin, pour N contacts d'amplitude inconnue, les premiers résultats sont encourageants. Actuellement, des cas plus complexes sont étudiés : superposition de plusieurs ondes avec des paramètres inconnus supplémentaires (nombres d'onde).

Figure 1 -

 1 Figure 1 -Modélisation du rail sur son support

Figure 2 -Figure 3 -

 23 Figure 2 -Atténuation des ondes de flexion verticales dans le rail

Figure 4 -Figure 5 -

 45 Figure 4 -Phase du champ acoustique rayonné par le rail à 350 Hz

ZcFigure 6 -Figure 7 -

 67 Figure 6 -Erreur d'estimation sur α en fonction de l'erreur sur l'atténuation, à 1000 Hz (∆ = 0.28 dB/m)

Figure 8 -

 8 Figure 8 -Visualisation de Λ autour du point solution (f = 550 Hz, z c = 3 m)

4. 4 . 2 Figure 9 -

 429 Figure 9 -Visualisation du Critère Λ (f = 550 Hz, z c1 = 3 m, z c2 = 5 m) Le critère Λ tel qu'il est défini présente un minimum plus marqué que lorsqu'il est défini en n'utilisant qu'une seule position pour l'antenne. Le point solution se situe à la croisée des vallées minimum de Λ 1 et Λ 2 . Les premiers tests effectués avec un algorithme d'optimisation basés sur une méthode SQP montrent que la convergence vers la solution est systématique ; le problème ainsi posé est mieux conditionné. De nombreux paramètres d'entrée peuvent et doivent être ajustés afin d'améliorer les performances de l'algorithme ; le choix des conditions initiales est notamment déterminant. L'objectif de notre étude ne portant pas sur cet aspect, nous nous arrêterons à ces résultats préliminaires obtenus pour une utilisation naïve de l'algorithme. Nous avons néanmoins pu vérifier qu'avec le critère (11), la convergence est assurée sur l'ensemble de la gamme de fréquence, même lorsque l'on initialise loin de la solution, ou proche d'une vallée minimum.

Table 1 -

 1 Elles sont caractérisées par une amplitude et un nombre d'onde complexe qui dépendent des paramètres physiques du rail et de son support (les paramètres utilisés dans cette étude sont donnés table1). On peut ainsi distinguer : -une onde potentiellement propagative d'amplitude F p et de nombre d'onde complexe k p ; -une onde de champ proche d'amplitude F d et de nombre d'onde complexe k d . Cette onde fortement atténuée est également appelée onde évanescente, elle existe dans un proche voisinage du point d'excitation. Paramètres mécaniques de la voie[START_REF] Hamet | Railway noise : use of the Timoshenko model in rail vibration studies[END_REF]

	Rail		Support
	Module d'élasticité (N/m 2 )	2 • 10 11	Masse d'une demie traverse (kg)	80
	Coefficient pertes internes	de	4 • 10 -3	Raideur de la se-melle (N/m)	3 • 10 8
	Module de cisaille-ment (N/m 2 )	7.7 • 10 10 Facteur de pertes dans la semelle	3.10 7
	Masse volumique (kg/m 3 )	8000	Raideur du ballast (N/m)	7.5 • 10 7
	Rigidité en flexion verticale (Nm 2 )	6.4 • 10 6	Facteur de pertes dans le ballast	3 • 10 7
	Masse par unité de longueur (kg/m)	60
	Coefficient de ci-saillement	0.4

1, on distingue deux types d'ondes de flexion qui peuvent se propager dans le rail, de part et d'autre de l'excitation.

Une source est placée en z 0 , puis tous les D zs mètres de part et d'autre du point d'excitation jusqu'à atteindre 60 dB d'atténuation pour l'onde vibratoire propagative. L'amplitude de chaque source est proportionnelle à la vitesse vibratoire du rail, elle-même déterminée en utilisant le modèle décrit dans la section 2.1. Ainsi, les sources élémentaires sont cohérentes et les champs associés interfèrent ; ceci confère au champ acoustique résultant des propriétés particulières (cf.[START_REF] Thompson | Investigation into the validity of two-dimensional models for sound radiation from waves in rails[END_REF]).

Méthode d'optimisation paramétriqueLe champ de pression décrit précédemment est mesuré au moyen d'un réseau de N mic microphones. Pour une fréquence f , on construit la matrice interspectrale Γ de dimension N mic × N mic . Parallèlement, une matrice spectrale modélisée Γ mod est calculée à partir du modèle vibro-acoustique décrit dans les sections 2.1 et 2.2. Elle dépend de plusieurs paramètres (amplitudes des forces de contact, nombres d'onde complexes dans le rail) que l'on cherche à déterminer par minimisation d'un critère basé sur l'erreur quadratique entre les matrices spectrales Γ et Γ mod :


Ruth Bergel
email: ruth.bergel@inrets.fr

H Rattaire
email: hugues.rattaire@free.fr

Maurice Aron
email: maurice.aron@inrets.fr

D Doucet
email: david.doucet@developpement-durable.gouv.fr

Eric Violette
email: eric.violette@developpement-durable.gouv.fr

Adeed risk in case of rain: some recent

Keywords: time series analysis, road safety, injury accidents, accident risk, added risk, rain

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

Introduction

Rainfall reduces the driver's vision and the friction of the vehicle on the roadway. The driver does not lower the speed of the vehicle accordingly, which induces a risk of running-off the road and a risk of rear end collision.

The researches led on these topics aim to improve the knowledge of the level of the additional accident risk in the case of rain, and, consequently, to implement safety preventive measures.

This paper is dedicated to the first objective -to improve the knowledge of the level of the added risk in case of rain. It provides aggregated results for the whole of France [START_REF] Rattaire | Sur-risque pluie en France métropolitaine[END_REF] which are compared to previous results obtained for a French region, the Haute-Normandie region, in the framework of the French project PREDITT (2004PREDITT ( -2008) ) IRCAD-SARI1 (Aron et al, 2006[START_REF] Bergel-Hayat | Added risk by rainy weather on the roads of Normandie-Centre region in France. Some elements for the French IRCAD-SARI traffic monitoring and information system for drivers and operators WCTR[END_REF], Following this introduction, a state of the art of the knowledge on the topic is briefly presented in section 2. The method for computing the injury accident risk and added-risk in case of rain is described in section 3. The estimated values for the accident risk and added risk are given in section 4, first for the Haute Normandie region and second for the whole of France. In section 5, some bias and limitations are discussed. We finally recall the results obtained, and give some research directions in order to improve them.

Previous results

After reviewing 14 studies addressing the relative risk in case of rain or snow, [START_REF] Violette | Weather as a Chronic Hazard for Road Transportation in Canadian Cities[END_REF] found a 70% increase in the risk of injury in case of rainfall in Canada. [START_REF] Keay | Road accidents and rainfall in a large Australian city[END_REF] noticed that, in Melbourne, the rain has an effect on the accident frequency and that this effect varies with the period in the day. They also demonstrated, as [START_REF] Lee | Analysis of Crash Precursors on Instrumented Freeways[END_REF] did before, that the amount of rainfall has a significant impact on accident occurrence.

In France, according to the French National Road Safety Observatory (ONISR, 2005), the proportion of road collisions during inclement weather was around 16% for the period 1995-2004, and 86% of these accidents occurred during rainfall. Rainfall is one of the weather variables used for monitoring the road safety level in the short term, both on a daily and monthly basis with the Giboulée model [START_REF] Bergel | A model of seasonally corrected indicators of road safety[END_REF] and is also explicitely used for explaining the road safety level in the medium term [START_REF] Bergel | Climate, road traffic and road risk -an aggregate approach[END_REF].

A review of time series modelling of the road safety level at national level confirms that, among all weather variables, rainfall is the most significantly correlated to the number of injury accidents and victims(.see for instance [START_REF] Scott | Modelling Time-Series of British Road Accident Data[END_REF], (Gaudry, Lassarre (Eds.), 2000), (Stipdonk(Ed.), 2008).

Researches have therefore focused on the value of his additional risk of injury accident (or added risk) due to rain. Previous studies stated that the added-risk in case of rain is about 2 ( [START_REF] Breton | Study of the rain impact on road safety (main roads and motorways). Statistical analysis[END_REF] or comprised between 2 and 3 [START_REF] Brodsky | Risk of road accident in rainy weather[END_REF]. In the case of France, studies led by both CETE of Normandie-Centre (the Technical Study and Engineering Centre of the French Ministry of Transport) and LCPC stated that the speed behaviour of drivers does not change sufficiently during adverse weather conditions, regarding the decrease of the braking performances of the vehicles (CETE, 2001). The availability of accurate and various weather data on the one hand, and of a huge traffic volume data set collected by CETE on the other hand, made it possible to estimate the value of the added-risk in case of rain accurately.

Method

3.1.

Risk and added risk formulas For a set of road sections A, and for a set of hours hH, let: Lr be the length of the road section r.

The risk 2 (number of injury accidents by vehicle kilometre) for that weather condition w is:

H h r r A h H h h r r r A h N w R w L w 1 1 ( ) ( ) ( ).           r Q (1)
The added risk in case of rain, noted ρ, is the ratio R (rain)/R (no-rain):

2 The term risk (in the place of "risk of injury accident") is used in this paper for commodity reasons.

Approximation

In (Saint-Pierre et al, 2007) a calculation is developed in order to take into account in K an incomplete information on the traffic variation in the case of rain., In this paper, the calculation will be done with the assumptions :

-either that there is no traffic variation in the case of rain. Then the ratio K is equal to the ratio duration norain

duration rain D _ _ 
. These durations are known in all places where meteorological stations are installed, and its is widely assumed that, due to the number of installed meteorological stations, rainfall duration is well known, apart from very particular places;

-or that the traffic decreases by δ which amounts to 1% or 2% during rain (figures validated on rural roads for light vehicles and trucks (Keays & Simmonds, 2005) (Khattak & de Palma, 1997). Thus K=(1-δ) D and an approximation of the added risk is given by:

D A ) 1 (     (3)
With these approximations, the calculation of the added risk is even possible when traffic data by rainy conditions are not available. Formulas (1), ( 2) and ( 3) can be applied to categories of accidents, categories of vehicles (if the traffic demand or its changes during rain periods is known by type of vehicle) and categories of networks. Even if the total length of elements of the road (bends, grades, intersections) is unknown and to compute the risk on these elements is thus not possible, the added risk by type of element is given by formula (3) -there are no more and no less curves on the road during rainfall than during no rain!

Results

The Haute-Normandie region

The 6 minutes rain gauge information from two French stations (Rouen and Evreux) has been validated with 1 minute human information f. A logistic regression model has been calibrated for each station, which led to increase the rain duration at Rouen and at Evreux. These adjusted meteorological data and the injury accident occurrence extracted from the French accident data files "BAAC" (Bulletin d'analyse des accidents corporels) have then been processed for the years 1995-2004.

By using hourly traffic counts available for the years 1997-2003 on 54 road sections of the region, the injury accident risk per vehicle-kilometre was estimated at 10 accidents for 100 millions of vehicle-kilometres during no rain periods, and at 20 accidents 100 millions of vehicle-kilometres during rain: the added risk, computed on these 54 counted sections, is thus close to 2 (Saint-Pierre &all, 2007). However, when computing it all road stations of the region (and not only on the 54 counted sections), the added-risk is estimated at.1.63 only; this decrease cannot be due only to the traffic approximation3 A good hypothesis, which will be confirmed from the results of the following section addressing the whole of France, is that traffic counts are not installed on very local roads where the risk is already high during no rain, and where the added risk4 is relatively low.

On bends, the added risk is higher (1,9 in the place of 1,63); among these accidents on bends, the added risk takes its highest value for front collisions (2.8), then for rear end or side-way collisions (2.2), then for isolated vehicles (1.6).

In intersection, the added risk is lower (1.5) However the preceding remark about the influence of the non-rainy and non-wet periods on the added risk value is relevant in the case of bends, since the water evacuation and the drying of the road must be shorter because of the elevation in curves, which leads to increase the value of the added risk in the case of rain. This remark also applies to the added risk value on grades.

The whole of France 5

The 6 minutes rain gauge information from about a hundred of French meteorological stations and the injury accident information6 extracted from the French accident data files "BAAC" have been processed for the years 1995-2008. Estimations of the added risk values for the whole of France were provided on a monthly basis over the period, and the raw values were still smoothed using a 12 months centered moving average filter (see Figure 1: Added risk in the case of rain for the whole of France for 1995-2008, raw data and smoothed data).

.A disaggregation according to several parameters (type of network, built-in/non built in areas, elements of the road) was also provided.

The whole of France: average

The added risk appears high (2.4 in average in 2004) 7 . The trend is a decrease, from 2, 8 in 1995, to 2, 4 en 2000 and to 2, 2 since 2005. This decrease is concomitant with a general decrease in the number of injury accidents, which has different causes: the road surface has been improved, the stability of vehicles has been enhanced, among others. The introduction of the speed enforcement policy, launched in France at the end of 2002, resulted in a significant decrease of the practised speed (see ONISR, 2005), and undoubtedly played an important role in these decreasing trends, particularly from 2003 onwards. Although the rain duration was especially high during the year 2003, the trend of the number of accidents in case of rain decreased and the added risk in case of rain significantly decreased [START_REF] Rattaire | Sur-risque pluie en France métropolitaine[END_REF] during the year.

4.2.2.

The whole of France: type of network In France as in other countries, motorways and local minor roads are at the two extrema of the security level (see Table 1: Risk by rain/no rain, according to several networks, for the years 2000 and 2005).

Table 1 provides the levels of risks for two particular years, before and after the speed enforcement implementation which took place in France in 2002. The decrease of the risk over time is clear. Considering the type of network, the level of the risk increases from motorways to main roads, and to local roads. As far as added risk in case of rain is concerned (see Table 2.Added risk according to the year and to the type of network, this range is not valid as there are two risk playing a role: the network and the rain. Indeed, in 2005 for instance, the added risk is the highest on main roads (2.64), second on secondary roads (2.41 on our "routes départementales") and on motorways (2.38), third on minor roads (1.88 on our "voies communales"). The added risk decreases from 2000 to 2005 by about 0.2 for every type of network, with the only exception of the secondary roads, for which the added risk value was unchanged. The whole of France: inside/outside built-in areas The added risk is higher outside built-in areas (2.6 in 2000 and since 2005) than the average value for the whole of France, and thus than inside built-in area (2 in 2005) (see.Table 3. Added risk according to the year and to inside/outside built-in areas.

It would be interesting to compute these two types of added risks according to the type of accident, to the period in the day, to the type of vehicles, in order to find an explanation to the differences of levels. Besides, the traffic flow is not very well known in built-in areas; and the assumption of an equal traffic whatever the meteorological condition need to be validated. . The decrease in the trend of the added risk in built-in areas (from 3 in 1995 to 2 in 2005) is much more significant than the decrease in the trend outside built-in areas (from 3 in 1995 to 2.6 in 2000 and 2005). The reason of such differences in the trends may be the respect of the speed limitation in urban areas.

4.2.4.

The whole of France: elements of the road It is well known that the added risk in case of rain varies according to some elements of the roadway (Shankar et all, 1995) (Shankar et all (2004). In France, the added risk reaches its highest value on bends (3.5 from 2005 onwards). On straight sections (outside intersection) its value is below the average (1.9 from 2005 onwards).

The decrease of the added risk between 1995 and 2005 is significant and about the same (0.9) on bends and on straight sections (outside intersection). The added risk value at intersections is intermediate (2.1 from 2005)(see Table 4.Added risk according to the year and to some elements of the road).

Discussion

The preceding results obtained for the whole of France are subject to bias as the rain duration is often under-estimated when it is measured with rain gauges. In addition, the influence of a wet roadway and of the presence of two-wheelers in the traffic flows need in all cases to be studied apart.

Nevertheless, the comparison of the added risk values according to several parameters, and the comparison of the trends of the added risk values for these categories are valid.

Possible bias

A rain gauge containing 0.2mm of water tips over each time it is full. However the rain duration may be under-estimated due to the evaporation, or to some block off of the gauge. A validation of a French rain gauge was carried out, and showed an important under-estimation.

With an under-estimation of c= %, the observed rain duration must be divided by (1-c). Thus the ratio D, defined in section 3.1, becomes:

) _ ( ) _ )( 1 ( _ duration rain c duration norain c duration rain     (4)
On the contrary, in the case of the Haute-Normandie region (see section 4.1), the computation of the added risk is not subject to this bias, as the measure of the rainfall was validated by human information.

Influence of a wet roadway

It is widely accepted that risk increases during no-rain periods for which the roadway surface is still wet. A number of roadway surface treatments, improving the friction, have been proved to reduce the risk of accident during no rain period for which the roadway surface is still wet. It should be reasonable, rather than matching the accidents during rain and no-rain periods, matching accidents during rain versus accidents during {no-rain and non-wet} periods in order to identify the effect of rain. Let us take an example. If the added risk during these periods were quite also important that during rain, and if the total duration in the year on the no-rain & non-wet periods were very little, the added risk in case of rain would systematically appear close to the value one.However it is very difficult to quantify the drying time of the roadway and thus to estimate the risk or the added-risk in these periods.

Influence of two-wheelers traffic

It is not possible to accept the assumption that the two-wheelers traffic does not decrease or decrease by δ= 1% or 2% only during rain. This assumption leads to an underestimation of the added risk in case of rain for two-wheelers, and therefore of the general added risk in case of rain. Moreover, as the percentage of traffic and accidents of two-wheelers increase in France, an apparent decrease of the added risk over time could be due to this assumption.

Conclusion and perspectives

The results confirm that the level of the added risk in case of rain is around 2, but may differ from 2 in a relevant manner.

The most detailed computation carried out on all rural roads of Haute-Normandie leads to a lower value (1, 6 only, whereas it is 2 when the computation is carried out on counted road sections only).

A more approximate computation carried out for the whole of France provides higher values (2.4 in average in 2004, the estimations ranging from 1.9 to 2.6 according to the network and road characteristics). These aggregate results underestimate the rain duration and thus over estimate the added risk in case of rain.

However, the analysis of the changes in the added risk value wrt time is valid.

A downward trend (the added risk decreased from 2.8 in 1995 to 2.2 from 2005 onwards) appears over the period 1995-2008. That relevant reduction in the level of the added risk in case of rain is concomitant with an increase in the security level, as the risk of accident decreased over the period.

Three research directions for improving these results can be mentioned. The influence of a wet roadway and of the two-wheeler' risk has to be taken account of. At last, the impact of rain on the injury accident severity remains to be studied.

7.

  volume for the hour h on the section r Q r be the total traffic volume on the section for the period H h N r (w) be the number of injury accidents for the hour h and weather condition w  be the proportion of the time during the hour h with the weather condition w h r (w)

Figure 1

 1 Figure 1 Added risk in the case of rain in France for 1995-2008 (raw and rmoothed data)

Table 1 . Risk (Number of Accident by 100 million vehicle kilometre) by rain/no rain, according to several networks, for the years 2000 and 2005

 1 More than 5,000 inhabitants. On built-in areas these figures must be taken cautiously, since the traffic levels on built-in areas are uncertain Source: BAAC, ONISR, Union Routière de France

	.
		Motorways	Main Roads Local Roads Built-in		Total
								Areas (*)
	Year	2000 2005 2000 2005 2000 2005 2000 2005 2000 2005
	Overall Risk	5	3	9	5	12	7	20	13	14	09
	No rain	5	3	8	4	11	6	18	12	13	08
	Rain	12	7	24	12	25	14	41	23	30	17
	(*)

Table 2 .Added risk according to the year and to the type of network 4.2.3.

 2

		Years End of 1995	2000	2005
	Autoroutes (Motorways)	3.5	2.58	2.38
	Routes Nationales	3.4	2.80	2.64
	(Main roads)
	Routes	départementales	2.9	2.42	2.41
	(Secondary roads)
	Voies Communales	2.5	2.09	1.88
	(Minor roads)
	Source: BAAC, ONISR.

Table 3 . Added risk according to the year and to inside/outside built-in areas.

 3

	Years End of 1995	2000	2005
	Added Risk
	France	3.0	2.38	2.18
	Built-in areas	2.8	2.24	1.95
	Outside Built-in areas(*)	3.2	2.56	2.55
	(*) outside built-in areas of more than 5,000 inhabitants
	Source: BAAC, ONISR.

Table 4 . Added risk according to the year and to some elements of the road

 4

	Years End of 1995	2000	2005
	Bends (outside inters.)	4.3	3.75	3.47
	Straight section (outside	2.8	2.05	1.86
	inters.)
	Straight section (inters.)	2.5	2.26	2.12
	Source: BAAC, ONISR.

The

[2004][2005][2006][2007][2008] French program SARI, funded by the French Ministry of Transport and managed by LCPC (the French Public Works Research Laboratory) aimed to install on the road network a warning sign system to alert drivers at moments there is a risk of running off the road, due to a poor adherence between the road and the tires. The sub-project IRCAD addressed the meteorological risk. Researches were led in three directions prior to the experiment: the drying of the roadway after the rain, the risk of running off the road according to the vehicle and driver's characteristics and according to the dynamic of the vehicle during rain, wind and wet roadway, and the added-risk in case of rain -seeAron et al., (2006[START_REF] Bergel-Hayat | Added risk by rainy weather on the roads of Normandie-Centre region in France. Some elements for the French IRCAD-SARI traffic monitoring and information system for drivers and operators WCTR[END_REF]. This third topic is discussed in this paper.

According to our approximation, there is no traffic decrease due to rain

According to our definition of the added risk

We consider the "Metropolitan" France, i.e. the historical France inside its natural borders.

Injury accidents with motorcycles have been disregarded here, since the decrease of motorcycles traffic during rain is likely to be significant but is not known

Acknowledgements

This work has been achievd in the followings of the sub-project "IRCAD" of the 2004-2008 French program "SARI", funded by the PREDIT and DRAST Departments of the French Ministry of transportation. We would like to thank Marie-Linne Gallene and Minh-Tan Do from LCPC who initiated and managed this project. The overall results would not have been obtained without the access to the weather-related information provided by to Meteo-France. Special thanks go to our colleague Jean-François Peytavin for his collaboration in collecting the accident-related information from the BAAC.

8.