Amelia Carolina Sparavigna

GIMP Retinex and Underwater Imaging

Keywords: Underwater Researches, Image processing, Retinex Filters. ______________________________________________

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

In some previous papers [START_REF] Sparavigna | Retinex filtering and thresholding of foggy images[END_REF][START_REF] Sparavigna | Effects of GIMP Retinex Filtering Evaluated by the Image Entropy[END_REF][START_REF] Marazzato | Retinex filtering of foggy images: generation of a bulk set with selection and ranking[END_REF][START_REF] Sparavigna | Evaluation of GIMP Retinex Filtering of Images by Means of the Shen++Max Shannon Entropy Finder[END_REF], we have discussed the use of the Retinex filter of GIMP, the GNU image manipulation program for X window systems, for enhancing the images recorded in foggy conditions and in night scenes. Since the GIMP Retinex filter is able to increase the visibility of objects in the abovementioned cases, such as in other circumstances [START_REF] Sparavigna | GIMP Retinex for enhancing images from microscopes[END_REF][START_REF] Sparavigna | An image processing approach based on Gnu Image Manipulation Program GIMP to the panoramic radiography[END_REF][START_REF] Sparavigna | The GIMP Retinex Filter Applied to the Fabric Fault Detection[END_REF], it is quite possible that it can be useful also for the underwater imaging. Here we will show some examples of the use of GIMP Retinex applied to images recorded in underwater conditions. As stresses in [START_REF] Zhang | Underwater image enhancement via extended multi-scale Retinex[END_REF], underwater exploration has become more active in recent years, with increasing applications, to name a few, in the studies of the marine life, for underwater engineering, search and rescue, and for the mapping of the seabed. In underwater exploration we have a vision that, due to the scattering and absorption effects which are giving hazing and bluish effects, results with low contrast, blur, and color variations [START_REF] Zhang | Underwater image enhancement via extended multi-scale Retinex[END_REF][START_REF] Sahu | A survey on underwater image enhancement techniques[END_REF]. On these premises, it seems that a Retinex filtering can help us for enhancing the quality of underwater images, as it is able to enhance the images of foggy and night scenes. Actually, several references exist which are discussing some Retinex algorithms for such a purpose [START_REF] Fu | A retinex-based enhancing approach for single underwater image[END_REF][START_REF] Zhang | Multi-scale retinex enhancement algorithm on luminance channel of color underwater image[END_REF][START_REF] Joshi | Quantification of retinex in enhancement of weather degraded images[END_REF][START_REF] Mirza | Enhancement of the Underwater Images Using Modified Retinex Algorithm[END_REF][START_REF] Raj | Underwater image enhancement using single scale retinex on a reconfigurable hardware[END_REF]; other methods for enhancing the underwater images are given in Ref. [START_REF] Iqbal | Underwater Image Enhancement Using An Integrated Colour Model[END_REF][START_REF] Wang | Saliency-based adaptive object extraction for color underwater images[END_REF][START_REF] Xinnan | Underwater Image Enhancement Algorithm Based on Contourlet Transform and Multi-Scale Retinex[END_REF][START_REF] Vasilescu | Color-accurate underwater imaging using perceptual adaptive illumination[END_REF][START_REF] Chiang | Underwater image enhancement by wavelength compensation and dehazing[END_REF][START_REF] Schettini | Underwater image processing: state of the art of restoration and image enhancement methods[END_REF][START_REF] Hitam | Mixture contrast limited adaptive histogram equalization for underwater image enhancement[END_REF][START_REF] Yang | Low complexity underwater image enhancement based on dark channel prior[END_REF][START_REF] Singh | Analysis of contrast enhancement techniques for underwater image[END_REF][START_REF] Hou | Automated underwater image restoration and retrieval of related optical properties[END_REF][START_REF] Kumar Rai | Underwater image segmentation using CLAHE enhancement and thresholding[END_REF][START_REF] Petit | Underwater image enhancement by attenuation inversion with quaternions[END_REF][START_REF] Yussof | Performing contrast limited adaptive histogram equalization technique on combined color models for underwater image enhancement[END_REF][START_REF] Andono | Underwater image enhancement using adaptive filtering for enhanced SIFT-based image matching[END_REF][START_REF] Ancuti | Enhancing underwater images and videos by fusion[END_REF][START_REF] Serikawa | Underwater image dehazing using joint trilateral filter[END_REF][START_REF] Prabhakar | An image based technique for enhancement of underwater images[END_REF][START_REF] Çelebi | Visual enhancement of underwater images using Empirical Mode Decomposition[END_REF][START_REF] Padmavathi | Comparison of filters used for underwater image pre-processing[END_REF][START_REF] Trucco | Self-tuning underwater image restoration[END_REF]. Before the discussion of the GIMP Retinex applied to underwater imaging, let us shortly remember the features of Retinex tools.

Retinex filters

Of the Retinex filtering, let us shortly remember the following. Retinex methods for image processing had been developed to solve some experimental observations. These experiments concerned the fact that it is easy to find discrepancies between an image that we have recorded by a camera and the real scene that we have observed. The reason is the following: humans are able to see details both in the shadows and in the nearby illuminated areas, whereas a photograph of the same scene is showing either the shadows as too dark or the bright areas as overexposed [START_REF] Barnard | Investigations into multi-scale Retinex[END_REF]. Modelled on the peculiar features of the human vision, which are quite different from those of the recording devices [START_REF] Barnard | Investigations into multi-scale Retinex[END_REF][START_REF] Bian | Retinex image enhancement techniques: Algorithm, application and advantages[END_REF][START_REF] Jobson | A Multi-Scale Retinex for bridging the gap between colour images and the human observation of scenes[END_REF], several algorithms of image processing had been developed to simulate the human ability of adapting to adverse conditions. These algorithms are known as Retinex algorithms. The first of them was conceived by Edwin H. Land, an American scientist and inventor, best known as cofounder of Polaroid Corporation. As explained in Ref. [START_REF] Barnard | Investigations into multi-scale Retinex[END_REF], through the years, Land evolved several models, until the last that he proposed in 1986. The term "Retinex" was coined by Land himself, combining the words "retina" and "cortex", to indicate the results of his researches. These researches tell that the human colour perception is involving all levels of vision processes, from the retina to the cerebral cortex. Several Retinex approaches exist [START_REF] Land | An alternative technique for the computation of the designator in the Retinex theory of color vision[END_REF][START_REF] Land | Recent advances in Retinex theory and some implications for cortical computations[END_REF][START_REF] Land | Recent advances in Retinex theory[END_REF][START_REF] Land | Experiments in color vision[END_REF][START_REF] Land | Color vision and the natural image[END_REF][START_REF] Jobson | Properties and performance of a center/surround Retinex[END_REF]: the single-scale Retinex (SSR), the multiscale Retinex (MSR), and, for colour images, the MultiScale Retinex with Colour Restoration (MSRCR). Among MSRCR we find the GIMP Retinex, a freely available tool developed by Fabien Pelisson [START_REF] Pelisson | GIMP Retinex[END_REF]. The resulting image of this filter can be adjusted selecting different levels, scales and dynamics. There are three "levels". The "uniform" level tends to process both low and high intensity areas fairly; the "low" level "flares up" the lower intensity areas on the image; the "high" level tends to "bury" the lower intensity areas in favor of a better rendering of the clearer areas of the image. The "scale" determines the depth of the Retinex scale. Minimum value is 16, a value providing gross, unrefined filtering. Maximum value is 250. Optimal and default value is 240. A "scale division" determines the number of iterations in the multiscale Retinex filter. The minimum required and recommended value is three. The "dynamic" slider allows adjusting colour saturation contamination around the new average colour (default value is 1,2).

Examples of the use of GIMP Retinex Images concerning scuba diving

Let us consider an underwater scene, depicted in the Figure 1. The original image is compared to the filtered images, with Uniform (U), High (H) and L (low) levels. In the Figure we see two different choices of the parameters. Let us note that, in this case, an overall effect of the filtering exists: it is that of enhancing the red tones. In the case of the example proposed in the Figure 2, the red tones are only slightly accentuated. Let us note that the use of the Retinex filter allows to see the details of the scuba divers and of the chain. Of course, other approaches for enhancing these images are possible. In the Figure 3 for instance, we used two other filters from GIMP on the original image of the Figure 2: by visual comparison we can conclude that the Retinex filter is giving a better result, in particular for the rendering of the environment where scubas are operating. On the left we can see image we can obtain using the "Equalize" command in the "Colors" menu of GIMP, applied to the original image of the Figure 2. This command "automatically adjusts the brightness of colors across the active layer so that the histogram for the Value channel is as nearly flat as possible, that is, so that each possible brightness value appears at about the same number of pixels as every other value. … Sometimes "Equalize" works very well to enhance the contrast in an image, bringing out details which were hard to see before. Other times, the results look very bad", as told in https://docs.gimp.org/en/gimp-layer-equalize.html. On the right, we can see the effect of the GIMP Cartoon filter, which modifies the image so that it looks like a cartoon drawing. It is clear from the Figure 4, that, in some cases, we are working as in a foggy scene [START_REF] Sparavigna | Retinex filtering and thresholding of foggy images[END_REF][START_REF] Marazzato | Retinex filtering of foggy images: generation of a bulk set with selection and ranking[END_REF]. The turbidity of the water is acting like the fog, diffusing the light and reducing the contrast. On the same basis, the GIMP Retinex is able to enhance the visibility of the underwater scene.

Marine Life

We can apply the GIMP Retinex for the study of the underwater images of the marine life. In the Figure 5 we are showing the picture of a whale (humpback whale). The individuals of this species are identified by the varying patterns on the tail flukes. In fact, the humpback whales have patterns of black and white pigmentation and scars on the underside of their tails that are unique to each whale [START_REF] Bank | Identifying Individual Humpback Whales[END_REF]. Other marks are scars, warts and bumps, pigmentation patterns or albinism, deformities, algal films or the shape and position of external features such as the dorsal fin [START_REF] Katona | Identifying humpback whales using their natural markings[END_REF]. Therefore photographic catalogues of whales exist and are maintained according to their marks [47]. The Figure 5 is showing that the use of GIMP Retinex can help in the identification of the individual by enhancing the visibility of the bumps on the skin. Other examples of the use of GIMP Retinex for the imaging of marine life are given in the Figures 6 and7. We can see that it is possible to detect the presence and count the individuals in the background of the picture (Figure 7).

Shipwrecks

The underwater wreck photography is difficult because of the large size of the subject and the limited light. In this case, a "black and white photography" is often preferred, because it is emphasizing the forms and the contrast, giving the possibility for enhancing the details in the wreck. Here some examples are given of the use of GIMP Retinex for an original colour image and for the same converted into grey tones. The Figure 8 shows, also in this case, that the Retinex filter is increasing the presence of red tones. The same happens if we use the "Equalize" command in the "Colors" menu of GIMP, as shown by the Figure 9. Actually, in the exploration of a shipwreck, the colour tones seem being not so relevant, because it is the presence of structures and other details which are important to give information about the wreck. Therefore, let us consider the original image of the Figure 8 and convert it into a grey-level image. The GIMP Retinex is working also on the images converted into grey tones and the result of the filter applied to the grey-tone image is remarkable (see the Figure 10). In the Figures 11 and12 we have the histograms of the two images in the Figure 10. The histograms are quite different, and therefore also the entropies of the images are different, being greater that of the filtered image [START_REF] Sparavigna | Evaluation of GIMP Retinex Filtering of Images by Means of the Shen++Max Shannon Entropy Finder[END_REF]. Two other examples are proposed in the Figures 13 and14. Figure 13 shows the result of the GIMP Retinex applied to an image of a research paper published in the Frontiers in Marine Science. The research work concerned in-situ corrosion measurements of WWII shipwrecks in Chuuk Lagoon. In this place, more than 65 shipwrecks and 250 aircraft were sunk during two main bombing raids. As told by [START_REF] Macleod | In-situ Corrosion Measurements of WWII Shipwrecks in Chuuk Lagoon, Quantification of Decay Mechanisms and Rates of Deterioration[END_REF], "the bunker fuel kept inside on-board storage tanks does present a real conservation management crisis". For this reason, an evaluation of corrosion rates is fundamental for the preservation of the local environment.

Underwater engineering. The engineering services, which are constituting a part of the underwater engineering, are mainly based on the works of the engineer divers. These divers have the capacities and knowledge needed to plan and execute the underwater operations concerning the specific engineering area. An example of engineering area which is requiring the engineer divers, is that of the civil infrastructures below water needing inspection and maintenance. In the case of underwater engineering works, a Retinex filtering of images could be useful for a detailed surveying of the structures and of the site of operations. In particular, the filtering of the images could be involved in a preliminary survey made by automated or remotely controlled devices, such as the underwater drones. In fact, the examples proposed above, and the others given in Figs. [START_REF] Iqbal | Underwater Image Enhancement Using An Integrated Colour Model[END_REF] and 16 for images in grey tones, seem evidencing that GIMP Retinex could be used to enhance the images recorded in preliminary surveying operations, to help the engineers in the planning of the following operations of the divers.

Figure 1 :

 1 Figure 1: In the left image, we can see the original picture (top/left) and the results of the GIMP Retinex filtering for Uniform (U), High (H) and Low (L) levels, for the parameters having the default values (Scale = 240; Scale division = 3; Dynamic slider = 1,2). (Original image, Courtesy jpatokal user of the wikitravel.org/en/User for Wikipedia). In the right image, we see again the original image and the results of the GIMP Retinex filtering for the three levels, in the case of parameters having values (Scale = 240; Scale division = 7; Dynamic slider = 3,0).

Figure 2 :

 2 Figure 2: Original image on the left and the Low level filtered image on the right, for the parameters having the default values (Scale = 240; Scale division = 3; Dynamic slider = 1,2). (Original image from flickr, Courtesy gingio for Wikipedia).

Figure 3 :

 3 Figure3: On the left we can see image we can obtain using the "Equalize" command in the "Colors" menu of GIMP, applied to the original image of the Figure2. This command "automatically adjusts the brightness of colors across the active layer so that the histogram for the Value channel is as nearly flat as possible, that is, so that each possible brightness value appears at about the same number of pixels as every other value. … Sometimes "Equalize" works very well to enhance the contrast in an image, bringing out details which were hard to see before. Other times, the results look very bad", as told in https://docs.gimp.org/en/gimp-layer-equalize.html. On the right, we can see the effect of the GIMP Cartoon filter, which modifies the image so that it looks like a cartoon drawing.

Figure 4 :

 4 Figure 4: Original image on the left and the Low level filtered image, for the parameters having values (Scale = 240; Scale division = 6; Dynamic slider = 2,0). (Courtesy University of Athens for Wikipedia).

Figure 5 :

 5 Figure 5: Original image and a detail of it enhanced using the Low level filtered image, for the parameters having values (Scale = 240; Scale division = 8; Dynamic slider = 4,0). (Original image, Courtesy Kurzon and Taragui for Wikipedia). We can see the eye and the skin covered with warts and bumps.

Figure 6 :

 6 Figure 6: The crystal jelly. On the left, the original image. In the middle the Low level filtered image for parameters (Scale = 240; Scale division = 3; Dynamic slider = 0,0), and on the right, for parameters (Scale = 240; Scale division = 3; Dynamic slider = 2,5). (Original image Courtesy Mnolf for Wikipedia).

Figure 7 :

 7 Figure 7: Aurelia jellyfishes. On the left, the original image and, in the right, the Low level filtered image for parameters (Scale = 250; Scale division = 7; Dynamic slider = 3,0). (Original image Courtesy Alexander Vasenin for Wikipedia). On the right, the Low level filtered image is given for parameters (Scale = 250; Scale division = 7; Dynamic slider = 0,0).

Figure 8 :

 8 Figure 8: The wreck of the Severance. On the left, the original image and, on the right, the Low level filtered image for parameters (Scale = 250; Scale division = 8; Dynamic slider = 3,0). (Original image Courtesy Viv Hamilton for Wikipedia).

Figure 9 :

 9 Figure 9: As in the case of the Figure 3, we show the results of the use of the "Equalize" command in the "Colors" menu of GIMP (left) and of the GIMP Cartoon filter (right). In fact, the image on the right is clearly showing that the rigging of the wreck.

Figure 10 :

 10 Figure 10: On the left, the original image has been converted into grey tones. On the right, the Low level filtered image for parameters (Scale = 250; Scale division = 8; Dynamic slider = 2,5). (Original image Courtesy Viv Hamilton for Wikipedia).

Figure 11 :

 11 Figure 11: Histogram of the original image in grey tones, from 0 to 255 (Figure 10, left).

Figure 12 :

 12 Figure 12: Histogram of the filtered image (Figure 10, on the right).

Figure 13 :

 13 Figure 13: On the shipwreck of the Yubae Maru. On the left the original image, in the middle the grey tone image and on the right the Low level filtered image for parameters (Scale = 240; Scale division = 8; Dynamic slider = 1,8). (Original image Courtesy from Article: MacLeod, I. D. (2016) In-situ Corrosion Measurements of WWII Shipwrecks in Chuuk Lagoon, Quantification of Decay Mechanisms and Rates of Deterioration. Front. Mar. Sci. 3:38. DOI: 10.3389/fmars.2016.00038).

Figure 14 :

 14 Figure 14: A shipwreck in the Tallinn Bay. Up/left, the original image and, on the right, the Low level filtered image for parameters (Scale = 250; Scale division = 8; Dynamic slider = 3,0). The same for the image in grey tones (Original image Courtesy Tuukritööde OÜ for Wikipedia).

Figure 15 :

 15 Figure 15: Up/left, the original image and, on the right, the Low level filtered image for parameters (Scale = 240; Scale division = 6; Dynamic slider = 2,0). The same for the image in grey tones (Original image Courtesy TPSDAVE, https://pixabay.com/it/tuffatore-immersioni-subacquee-nuoto-85913/).

Figure 16 :

 16 Figure 16: Up/left, the original image and, on the right, the Low level filtered image for parameters (Scale = 240; Scale division = 6; Dynamic slider = 2,0). The same for the image in grey tones (Original image Courtesy TPSDAVE, https://pixabay.com/it/mare-ocean-acqua-subacquea-luce-86425/).


Alfonso M Panunzio

G Puel

R Cottereau
email: regis.cottereau@centralesupelec.fr

S Simon

X Quost

Construction of a stochastic model of track geometry irregularities and validation through experimental measurements of dynamic loading

Keywords: track geometry irregularities, railway vehicle dynamics, statistical validation

This paper describes the construction of a stochastic model of urban railway track geometry irregularities, based on experimental data. The considered irregularities are track gauge, superelevation, horizontal and vertical curvatures. They are modelled as random fields whose statistical properties are extracted from a large set of on-track measurements of the geometry of an urban railway network. About 300 to 1000 terms are used in the Karhunen-Loève/Polynomial Chaos expansions to represent the random fields with appropriate accuracy. The construction of the random fields is then validated by comparing on-track measurements of the contact forces and numerical dynamics simulations for different operational conditions (train velocity and car load) and horizontal layouts (alignment, curve). The dynamics simulations are performed both with and without randomly generated geometrical irregularities for the track. The power spectrum densities obtained from the dynamics simulations with the model of geometrical irregularities compare extremely well with those obtained from the experimental contact forces. Without irregularities, the spectrum is 10 to 50 dB too low.

Introduction

The layout of a railway network is decided in the design phase according to engineering considerations. However the track evolves during its operational life because of train loads, weather conditions, substructure variations and maintenance interventions. In particular, track geometry irregularities appear and impact the dynamical response of a train [START_REF] Lei | Analyses of dynamic response of vehicle and track coupling system with random irregularity of track vertical profile[END_REF][START_REF] Gullers | High-frequency vertical wheel-rail contact forces-Field measurements and influence of track irregularities[END_REF][START_REF] Zhai | Fundamentals of vehicle-track coupled dynamics[END_REF][START_REF] Perrin | Quantification of the influence of the track geometry variability on the train dynamics[END_REF][START_REF] Yang | Effect of track irregularity on the dynamic response of a slab track under a high-speed train based on the composite track element method[END_REF], the rolling contact fatigue [START_REF] Evans | Vehicle / track interaction and rolling contact fatigue in rails in the UK[END_REF][START_REF] Li | An investigation into the causes of squats Correlation analysis and numerical modeling[END_REF][START_REF] Karttunen | The influence of track geometry irregularities on rolling contact fatigue[END_REF] and the rail corrugation [START_REF] Jin | Effect of track irregularities on initiation and evolution of rail corrugation[END_REF]. In fact, even a small variation of the geometry can produce large changes of the rail-wheel contact patch and forces [START_REF] Sun | Effect of track geometry irregularities on wheel-rail impact forces[END_REF]. The maintenance policy of a railway line is very sensitive to track irregularities, especially on an urban line that is characterized by intense traffic and strongly variable operational conditions (velocity and passenger loads). Therefore, the track geometry is constantly monitored using different measurement systems [START_REF] Grassie | Measurement of railhead longitudinal profiles: a comparison of different techniques[END_REF][START_REF] Fröhling | Measurement, interpretation and classification of South African track geometry[END_REF][START_REF] Haigermoser | Road and track irregularities: measurement, assessment and simulation[END_REF][START_REF] Nielsen | Overview of Methods for Measurement of Track Irregularities Important for Ground-Borne Vibration[END_REF].

Using information from measurements, it is possible to estimate the statistical properties of the track geometry and generate the irregularities as random fields indexed on space. The first step consists in extracting realizations of the fields from the measurements set and estimate the probability density functions and the spatial dependency (usually represented by the power spectral densities, PSDs). In many railway applications the irregularities are generated from the PSD as Gaussian fields [START_REF] Karttunen | The influence of track geometry irregularities on rolling contact fatigue[END_REF][START_REF] Claus | Modeling and Simulation of Railway Bogie Structural Vibrations[END_REF][START_REF] Iyengar | Random field modeling of track irregularities[END_REF][START_REF] Zhu | Advanced Stochastic Modeling of Railway Track Irregularities[END_REF][START_REF] Majka | Dynamic response of bridges to moving trains: A study on effects of random track irregularities and bridge skewness[END_REF] using a spectral representation [START_REF] Shinozuka | Simulation of stochastic processes by spectral representation[END_REF]. In [START_REF] Perrin | Track irregularities stochastic modeling[END_REF] the geometry irregularities of a high-speed railway network have been generated as a non-Gaussian irregularities random field with a more complete dependence structure. The method proposed in [START_REF] Perrin | Track irregularities stochastic modeling[END_REF] combines the Karhunen-Loève Expansion [START_REF] Loève | Probability theory[END_REF], to capture the spatial dependency, and the Polynomial Chaos Expansion [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] to obtain a good statistical representation of the non-Gaussian fields. This technique is used in this work to generate track geometries which are statistically representative of the RATP urban railway network (RER line A), operating in the Paris metropolitan area. The main differences between a high-speed and an urban railway network are the operational conditions (velocity, passengers loads and train passages frequency) and horizontal layout (curve frequency and radius). This implies a difference in the interactions between vehicles and track and the maintenance policies, and so the geometric irregularities. The main objective of this paper is to construct a stochastic model of the track geometry irregularities using experimental measurements.

In a second step the stochastic model is validated employing on-track wheel-rail forces measurements. The measurements of the geometry irregularities can be directly used as input for deterministic numerical computations to simulate the passage of the train on a real track portion [START_REF] Pombo | An alternative method to include track irregularities in railway vehicle dynamic analyses[END_REF][START_REF] Luber | Track geometry evaluation method based on vehicle response prediction[END_REF][START_REF] Li | Assessing track geometry quality based on wavelength spectra and track-vehicle dynamic interaction[END_REF]. Dynamic simulations, using the constructed random irregularities generator, are performed for different operational conditions using the multibody simulation software SIMPACK. A statistic characterisation of the wheel-rail contact forces is therefore obtained. For the validation of the irregularities model, the results of dynamic simulations are statistically compared with measurements of dynamic quantities of interest. These measurements are performed on a vehicle equipped with strain gauge systems positioned in the wheels [START_REF] Gullers | High-frequency vertical wheel-rail contact forces-Field measurements and influence of track irregularities[END_REF][START_REF] Kanehara | Measuring rail/wheel contact points of running railway vehicles[END_REF]. A comparison of the power spectral densities of the measured and the simulated contact forces is made to validate the stochastic model of the track irregularities.

The paper is organized into two parts: the Section 2 is dedicated to the track geometry measurement and random modelling, while in Section 3, in order to validate the model, the spectra of the simulated and experimental dynamic wheel-rail forces are compared.

Stochastic model of the track irregularities

The measurements used to extract the statistical features are firstly described in this section. Then, the statistics of the track irregularities and the technique used to construct a stochastic model of them are presented.

Track geometry measurement data description

By considering the rails as lines evolving in the space, the track geometry is the set of fields defining their position. Different choices can be made on these fields, but it is common, in the railway domain, to work with:

• Track gauge, G(s), the spacing of the rails measured between the inner faces • Superelevation, E(s), the difference in height between the two rails • Horizontal curvature, C H (s), referred to the track median line • Vertical curvature, C V (s), referred to the track median line where the indexation variable s is the curvilinear abscissa of track median line. A geometry irregularity is the difference between the design and the current geometry. In some cases, instead of working with the curvatures irregularities, one can introduce two quantities: alignment (al(s)) and vertical profile (lp(s)). These quantities are defined as the double integral in space of the curvatures irregularities (C Hirr (s) and C V irr (s)) and represent the displacement of the mean track line from its design position:

al(s) = C Hirr (s)ds ; lp(s) = C V irr (s)ds (1)
The geometry irregularities are sketched in Figure 1.

Measurement system

The measurement train V355 system consists in a mobile rigid base equipped by two fibre optic gyroscopes and six laser cameras pointing the rails (three on each side). It is represented in Appendix A. The superelevation and the gauge are directly measured whereas the horizontal and vertical curvatures are obtained using a chord-based system [START_REF] Grassie | Measurement of railhead longitudinal profiles: a comparison of different techniques[END_REF]. This means that the curvatures are not directly measured, but they are linked to a versine measurement, v(s), with a relation depending on the spatial frequency, f , that can be modelled with a transfer function, H(f ), and its associated impulsive response, h(s) [START_REF] Aknin | A new approach for the modelling of track geometry recording vehicles and the deconvolution of versine measurements[END_REF][START_REF] Mauer | Determination of Track Irregularities and Stiffness Parameters with Inverse Transfer Functions of Track Recording Vehicles[END_REF]:

H(f ) = 1 (2πf ) 2 1 - A 2 A 1 + A 2 exp(-2A 1 πf i) - A 1 A 1 + A 2 exp(2A 2 πf i) v(s) = (h * C H )(s) h(s) = F -1 (H(f )) (2)
December 17, 2016 Vehicle System Dynamics paper where A 1 and A 2 are the distances between the laser cameras (4.5 m and 5.5 m) and the symbol * denotes the convolution product. The phase of this transfer function, shown in Figure 2, depends on the circulation direction because of the asymmetry of the system (A 1 = A 2 ). To avoid a noise amplification, the transfer function in Equation 2 is not brutally inverted but a Wiener deconvolution approach is employed [START_REF] Aknin | A new approach for the modelling of track geometry recording vehicles and the deconvolution of versine measurements[END_REF][START_REF] Panunzio | A stochastic model of the railway track geometry[END_REF].

Frequency [m -1 ] 0 0.5 1 1.5 2
Module [m], in dB 2): (left) module and (right) phase.

Measurement campaigns

All the measurement campaigns from 2005 to 2015 concerning the RER A railway network have been processed. This is an urban railway network crossing the Paris metropolitan area West to East (75 km double-track). Only the ballasted portions (93 % of the total length) are processed in this paper. The maximal operational velocity goes from 80 to 120 km/h according to the portion. With around 350 million of passengers per year, it is one of the busiest network in the world. During the considered period, the measurement train passed about 22 times on each point of the network. The passages of the train on the secondary (storage and intermediate) tracks and points and crossings have been suppressed. At the end about 9 millions of measurement points, recorded with a sampling rate equal to 0.25 m, are considered for the statistical model. An example of the track geometry measures with the corresponding theoretical (design) values is shown in Figure 3.

Curvilinear abscissa [m]

Statistical description of the track geometry measurements

After collecting the measurements data, realizations of the geometry fields can be extracted from the measurement campaigns. The irregularities of superelevation, horizontal and vertical curvatures are defined as the differences between the measurements and the corresponding design (or theoretical) known values:

     E(s) = E d (s) + E irr (s) C H (s) = C Hd (s) + C Hirr (s) C V (s) = C V d (s) + C V irr (s) (3)
where the subscripts d and irr indicate respectively the design values and the irregularities fields.

Concerning the gauge, the design value is equal to 1435 mm on the whole line. However the measurements show that the gauge tends to increase in the curved track portions. In fact, the correlation coefficient between the gauge measurements and the horizontal design curvature |C Hd (s)| (taken in absolute value to equally consider leftward and rightward curves) is equal to 0.57 (while for the other irregularities it is lower than 0.02). The correlation table of the irregularities fields and the design curvature is shown in Figure 4. Therefore the gauge field is separated into two terms. One of them, G 0 (|C Hd (s)|), depends on the design curvature and is calculated by linear regression, while the other one, G irr (s), is the irregularity field that is decorrelated from the curvature and has to be modelled:

G(s) = G d + G 0 (|C Hd (s)|) + G irr (s) G 0 (|C Hd (s)|) = 3.97|C Hd (s)| mm•km (4)
with the curvature expressed in km -1 . The irregularities fields thus defined are centred and decorrelated from the horizontal layout of the track. Their experimental histogram, here used as probability densities functions (PDFs), are shown in Figure 5. The marginal PDFs of the fields are not Gaussian.

The cross-correlations between the irregularities are small (lower than 0.05 at zero lag, Figure 4). The experimental joint PDFs (at zero lag), Figure 6, show that there are not strong dependencies between the fields. Therefore, the dependences between the irregularities are neglected in this work and the fields are identified separately.

The experimental correlation functions, Figure 7, show that the correlation length of each field is different: around 200 m for the superelevation and the gauge, less than 20 m for both the curvatures. The corresponding experimental power spectral densities (PSDs) are shown in Figure 8. The superelevation and the gauge correlation functions, which have a decreasing exponential behaviour, show that these irregularities are more related to the topography of the track. Similar irregularities PSDs can be found in [START_REF] Karttunen | The influence of track geometry irregularities on rolling contact fatigue[END_REF][START_REF] Zhu | Advanced Stochastic Modeling of Railway Track Irregularities[END_REF][START_REF] Liu | Characteristic Analysis of the PSD of Track Irregularity on Ballastless Track in High-speed Railway[END_REF]. Note however that the correlation length of the horizontal design track curvature is more than three time larger (Figure 9) and shows that there are not components of the irregularities correlated with the horizontal layout. The curvatures irregularities are more related to the rail short lengths such as corrugation and welds [START_REF] Iyengar | Random field modeling of track irregularities[END_REF][START_REF] Weston | Monitoring vertical track irregularity from in-service railway vehicles[END_REF][START_REF] Weston | Monitoring lateral track irregularity from in-service railway vehicles[END_REF]. The histogram of the length between two welds (on the considered railway network) superposed with the correlation of the horizontal curvature is shown in Figure 10. Note the correspondence between the histogram and the correlation function peaks. In the literature, the alignment al and the vertical profile lp (Equation 1) are often used instead of the horizontal and vertical curvatures. The PSDs of these quantities, obtained from the curvatures measurements, are shown in Figure 11 with the two German PSDs analytical standard that are usually used for simulation in European countries [START_REF] Zhang | Non-stationary random vibration analysis for train-bridge systems subjected to horizontal earthquakes[END_REF][START_REF] Lin | Theory analysis and test research of Chinese main track irregularities PSD[END_REF]. A strong correspondece between them can be observed especially in the length range from 3 m to 40 m. The weather influence on the track geometry has been analysed by estimating their statistics taking the measurements in the summer and winter seasons separately. No differences have been observed. That is why the weather is not considered into the model. Around 35 % of the line is in tunnels where no significant temperature variations are registered. The rest of the line is essentially composed by continuous welded rails, in which only the breathing zones can be thermally deformed along the longitudinal direction. This deformation does not influence the track geometry (rails positions in the lateral and vertical directions).

Lag [m]

Parametric model of the random fields

In this paper the fields generator is constructed by combining the Karhunen-Loève expansion (KLE) and the Polynomial Chaos Expansion (PCE). This strategy has been successfully applied to model the track irregularities of a high speed railway network [START_REF] Perrin | Track irregularities stochastic modeling[END_REF].

Karhunen-Loève expansion

Using a Karhunen-Loève expansion [START_REF] Loève | Probability theory[END_REF], a random field f (s) (that could be any irregularity field) is approximated by its truncated projection on an orthogonal basis. In this way it is expressed as a sum of M deterministic spatial (depending only on the curvilinear abscissa) functions ϕ i (s) multiplied by random decorrelated and normalized coefficients η i :

f (s) ≈ M i=1 λ i ϕ i (s)η i L 0 C(s, s )ϕ i (s )ds = λ i ϕ i (s) E[η i η k ] = δ ik (5)
where s, s ∈ [0, L], C(s, s ) is the covariance function associated with the field, δ ik is the Kronecker delta and E denotes the expectation operation. In a more compact form, this expansion is written as:

f (s) ≈ Φ(s) √ ΛH (6)
where Φ(s) is a M -dimensional row vector such that Φ l (s) = ϕ l (s), Λ is a M × M diagonal matrix such that Λ lm = δ lm λ l and H is a M -dimensional column vector such that H l = η l December 17, 2016 Vehicle System Dynamics paper This expansion reduces the normalized mean-square error ε KL resulting of its truncation that, thanks to the orthogonality property of the basis, can be expressed as:

ε KL = 1 - M i=1 λ i E[f (s) 2 ] = (S(f ) -S g (f ))df S(f )df (7)
where S(f ) and S g (f ) are respectively the reference (estimated from measurements) and the modelled PSDs. The error on the PSDs (chosen equal to 1%) only depends on the truncation error of the KLE. The length L in Equation 5is different for each irregularity field because of the different correlation lengths (Figure 7). It is set to three times the correlation length. Taking a larger length does not change the statistical properties but would require more terms in the KLE. The threshold for the truncation error is set to 1% in this work. This error is evaluated for the four fields separately (Figure 12). The field which needs more terms to be identified is the gauge (314 terms to obtain a truncation error equal to 0.01). The first (most energetic) three modes (deterministic functions ϕ i (s)) associated with each field are shown in Figure 13. Note the different lengths of the modes according to the field.

No. of terms [-]

Polynomial Chaos expansion

If a random field is Gaussian, then the coefficients η i are independent standard Gaussian variables. In a general case their distribution has to be identified. The multivariate distributions of the random projection coefficients (η i ) of the basis of each random field are characterized using the PCE [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]. This consists in expanding the random coefficients in a sum of deterministic coefficients β ij multiplied by a polynomial basis ψ j of multivariate

η i ≈ N j=1 β ij ψ j (ξ) (8)
The polynomial functions are chosen such that they are orthogonal with respect to the distributions of the germs. For instance, Hermite polynomials are chosen with normal germs while Legendre polynomials are chosen with uniform germs. Given the number of germs n g and the polynomial degree d the number of polynomials N is:

N = (n g + d)! n g !d! (9)
In a more compact form, the Equation 8 is written as:

H = BΨ (ξ) (10)
where

Ψ (ξ) is a N -dimensional column vector such that Ψ l (ξ) = ψ l (ξ) and B is a M ×N matrix such that B lm = β lm .
Since the number of coefficients (dimension of the random space to characterize) is very high, advanced techniques to identify the deterministic coefficients are employed to guarantee a good representation [START_REF] Soize | Identification of high-dimension polynomial chaos expansions with random coefficients for non-Gaussian tensor-valued random fields using partial and limited experimental data[END_REF][START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF]. These coefficients are thus chosen as the argument that maximizes the log-likelihood between the available realisations of the random coefficients and their chaos representation. Because of the decorrelation and the normalization properties due to the KLE (Equation 5) the following constraint has to be satisfied:

BB T = I ( 11
)
where I is the identity matrix.

For each KLE coefficient the Kullback-Leibler distance [START_REF] Kullback | On information and sufficiency[END_REF], evaluated in the range of the 0.1th l and the 99.9th u percentiles of the target PDF, can be used to define an error i between the target (obtained from measurements) and the generated PDFs. The maximal value of this error, on all the KLE coefficients, is taken as PCE identification error ε P CE :

i = u l p(x) ln p(x) q(x) dx ε P CE = max( i ) (12)
where p(x) and q(x) are the target and the generated PDFs. Choosing a threshold tolerance τ (0.01 in this work), the number of germs and the polynomial degree are chosen as:

(n g , d) = arg min(N (n g , d) | ε P CE < τ ) (13)
For instance, in Figure 14, the value of the error ε P CE defined in Equation 12 is shown for different numbers of germs and polynomial degrees in the case of the vertical curvature. Note that, when the number of PCE terms is lower than the number of KLE terms to model (N < M ), the identification is not possible since the Equation 11 cannot be satisfied. 12) of the vertical curvature irregularity (100 KLE terms). The dashed red line delimits the area in which the identification is not possible since the number of PCE terms is lower than the number of KLE terms to model. The dotted white line delimits the area in which the PCE error is lower than 1%.

The number of KLE and PCE terms (with Hermite and Legendre polynomial basis) for each field are shown in Table 1. Using Legendre polynomials (and so uniform germs) requires more terms to obtain the same accuracy. Therefore, Hermite polynomials (and so Gaussian germs) are chosen for the model.

Random generation

The germs distribution determines the polynomial family used as basis. In this work Gaussian germs and so Hermite polynomials are employed. Finally, by generating ν germs ξ it is possible to obtain ν independent realisations of the generated field f gen (s):

f gen (s) = f d (s) + f 0 (|C Hd (s)|) + M i=1 N j=1 λ i ϕ i (s)β ij ψ j (ξ) (14)
where f d (s) and f 0 (|C Hd (s)|) are respectively the design value and the component, not null for the gauge, correlated to the design horizontal curvature.

Generation with an arbitrary length

The length of the generated fields, L, is set by the integration interval on which the Karhunen-Loève expansion is defined (Equation 5). However, one usually needs to generate a field longer (L new ) than that length. The simple concatenation of more realisations will break the continuity of the field and the correlation structure around the breaking points. The generation method described in [START_REF] Cho | Karhunen-Loève expansion for multi-correlated stochastic processes[END_REF], originally proposed to generate crosscorrelated random fields, is here adapted to be used to outcome this issue in the case of stationary random fields. The main steps of the adapted method are here summarized:

• Define the function:

C(|s -t|) = C(|s -t|), if |s -t| ≤ L 0, if |s -t| > L (15)
where C(|s -t|) is the covariance function of the considered field. • Define the M × M matrix K such that:

K lm = 1 √ λ l λ m L 0 2L L C(|s -t|)ϕ l (s)ϕ m (t -L)dsdt = 1 √ λ l λ m L 0 L 0 C(|s -t -L|)ϕ l (s)ϕ m (t)dsdt (16)
• Generate, as in Equation 10, a n-dimensional set (such that nL > L new ) of KLE coefficients: H (k) = BΨ (ξ (k) ), with k going from 1 to n. • Define the 2nM × 2nM matrix R such that:

RR T =         I K • • • 0 K T I K . . . . . . . . . . . . . . . K T I K 0 • • • K T I         (17)
where 0 is a M × M null matrix. • Define a n-dimensional set of transformed KLE coefficients H such that:

        H(1) . . . H(k) . . . H(n)         = R        
H (1) . . .

H (k)

. . .

H (n)         (18)
• Finally, the field g(s) (of length L new ) is generated as:

g(s) =                  Φ(s) √ Λ H(1) , if 0 ≤ s ≤ L . . . . . . Φ(s -(k -1)L) √ Λ H(k) , if (k -1)L < s ≤ kL . . . . . . Φ(s -(n -1)L) √ Λ H(n) , if (n -1)L < s ≤ L new (19)
This method works under the condition that, for a spatial lag larger than L, the correlation function of the field can be approximated to zero. The covariance function of the field g(s) is equal to C(|s -t|) (defined in Equation 15) on all the domain. This is due to the fact that, by construction:

E[ H(k) ] = 0 E[ H(k) ( H(l) ) T ] =        I, if k = l K, if k = l -1 K T , if k = l + 1 0, otherwise (20)
where 0 is a M -dimensional null column vector.

Even if the field g(s) is defined as piecewise, its continuity is guaranteed also at the breaking points when M → ∞. For a discussion about the continuity see Appendix B.

Comparison with measurements

The KLE allows to capture the spatial correlation of the fields and the PCE ensures the statistical variability of the track geometry: the obtained realisations are representative of a whole railway network. The PSDs and the marginal distribution of the four fields are shown in Figure 15 and Figure 16. The frequencies that are not properly represented have an energetic level about 20 dB lower than the maximal power. With a KLE truncation error chosen equal to 1%, the curvature fields are correctly represented in the spectrum range of interest (wavelengths between 1 m and 150 m ), differently from the gauge and the superelevation. Although taking more terms would decrease the error, the effect of more terms on the dynamic forces spectra is limited (as shown in Section 3.3). The error on the distributions, that is evaluated using the Kullback-Leibler distance, depends on the PCE identification and is indicated in Figure 16. It is shown that the December 17, 2016 Vehicle System Dynamics paper model is able the capture the distribution tails. Note that the condition used in Equation 13 concerns the maximal distance of all the KLE coefficients and not the distance on the marginal PDFs of the field. This means that all the modes, even the less energetic ones, are properly represented. Some examples of measured and generated irregularities are shown in Figure 17.

Wavelength [m]

Gauge [mm]

Curvilinear abscissa [m]

Dynamic simulations and validation

In this section, track irregularities stochastic model, described in Section 2, is validated through dynamical simulations. Measurements of the wheel-rail contact forces are collected on the same railway network on which the random generator is built. The goal is to compare the experimental spectra of the contact forces with the ones obtained by simulating the passage of the train on tracks modelled with random irregularities.

Wheel-rail contact forces experimental data

Measurement system

The wheel-rail contact forces are usually decomposed in the track frame into vertical (Q), lateral (Y) and longitudinal (X) directions (Figure 18). An adapted train running on the RER rail network (MI 2N unit) is used to measure the wheel-rail contact forces. The leading passenger car and the motor coach, which is positioned behind the passenger car, are instrumented. These vehicles have two bogies and four wheelsets as sketched in Figure 19. The wheels of the front bogies are drilled and equipped with strain gauges (Figure 3.1.1) able to measure the deformations which are processed to obtain the contact forces. The rear bogies are equipped with strain gauges located on the axle boxes. The measurement of the forces on the rear bogies is less precise than the front ones, because the gauges are more distant from the wheel-rail contact interface and they are filtered by December 17, 2016 Vehicle System Dynamics paper the wheelset (that is not infinitely rigid). In this work the validation will be performed by using the measures from the strain gauges located in the wheels (front bogies). The angle between the two cars is measured and used to correctly locate the train on the track. In fact the angle signal is strongly correlated with the horizontal curvature of the track, i.e. the spatial lag between angle and horizontal curvature is given by the peak of their cross-correlation function. An example of the measured forces on a track portion with three curves is shown in Figure 21. As observed in this figure, the curve radius has a big impact on the forces fluctuations and their mean values (especially for the front wheelset which guides the bogie in the curve inscription) introducing a strong spatial non-stationarity. Examples of wheel-rail contact positions in a curved track are shown in Figure 22.

Measurement campaigns

The measurements are performed on the RER line A (the same network on which the geometry is measured by the other recording train) at different operating conditions (maximal running velocity going from 80 to 120 km/h), in order to be representative of the trains circulating the this railway network.Two levels of weight of the leading car are used in the campaigns to simulate the passenger load: empty and charged (6 passengers per m 2 , i.e. 20960 kg). The totality of the measurements corresponds to 1172 km of track travelled (630 km with an empty leading car and 542 with a charged one).

Comparison of the power spectral densities of the contact forces

The measurements of the dynamic forces are compared with numerical simulations where the track geometry are generated as in Section 2.4. To compute the dynamical forces, an advanced semi-Hertzian contact model [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF] is used in the simulations. In this model the rail profile is discretized in strips along the lateral direction. Then, in each strip, the contact problem is solved as Hertzian. The railway dynamics software SIMPACK is used for this purpose. The geometry irregularities are the only variable parameters in the simulations. In reality, other track parameters should be considered as stochastic to refine the simulations, but, since the objective is to validate the geometry irregularities model, these other parameters are used as deterministic. The rail profile (UIC60) is considered as theoretical, i.e. the wear is not included. The wheel-rail friction coefficient varies neither along the rail profile nor along the curvilinear abscissa of the track: it is equal to 0.2 in all the simulations. The rail cant is set to 1/20 (standard in France) and does not vary dynamically.

A multi-body model representative of the passenger car of the MI 2N train (circulating on the RER A railway network) is used for the simulations. The vehicle is composed by a car enclosure supported by two bogies equipped by two wheelsets each. The stiffness and the damping of the primary and secondary suspensions are set to the manufacturer data in all the simulations. The curb weight of this vehicle is approximatively equal to 55 t. To take into account the weight of the measurement system, 4 t uniformly distributed in the car enclosure have been added.

At this stage the objective is to validate the geometry irregularities model by comparing numerical simulations to measurements. For this purpose all the measurement campaigns (Section 3.1.2) are replicated by simulating the passage of the vehicle on a track affected by geometry irregularities. Therefore, in each simulation the track design (nominal curvatures, gauge and superelevation), the weight and the running velocity are the same as in the corresponding measurement campaign, while the geometry irregularities are randomly generated according to the model presented in Section 2.4.

In this section we analyse the spectral densities of the contact forces in the length range from 1 m to 150 m (largest correlation length of the irregularities fields). The train dynamics strongly changes according to the curve radius (and the associated superelevation) and the velocity, as shown in Figure 21. For this reason, the PSDs are separately computed in straight and curved (outer and inner rails) track portions for constant velocities, acceleration and breaking phases.

Straight track

The PSDs of the contact forces on the front wheelset in straight track portions rolling at constant velocity are shown in Figure 23. When the geometry irregularities are taken into account in the simulations, a strong improvement is produced on the PSDs for the three components of the contact force. In the simulations performed without irregularities, the spectra of the contact forces are between 1 and 2 orders of magnitude smaller.

Some peaks appear, especially on the PSD of the longitudinal forces X, at a length corresponding to the wheel circumference and its submultiples. The wheel nominal diameter is equal to 0.42 m. These peaks seem to be due to wheel untrueness (see Section 3.4).

Another discrepancy, observable on the spectrum of the longitudinal force X between the lengths of 7 m and 50 m, could be due to the train traction system which is not considered in the simulations. As shown in [START_REF] Simon | De la dynamique ferroviaire à l'accommodation microstructurale du rail Phd thesis[END_REF], the motor traction effects affect the longitudinal equilibrium of the train in alignment by modifying the contact location and so the rolling radius.

Curved track

The PSDs of the contact forces on the front wheelset in curved track portions rolling at constant velocity are shown in Figure 24. In a curved track the variability of the contact forces is larger than a straight track especially for the lateral Y forces on the outer rail.

In fact, due to the localisation of the contact on the wheel flange, a small variation of the wheel-rail position causes larger changes of the contact force. When the irregularities are not considered the contact position is more stable. The random model of the irregularities allows to very well represent the forces spectra on all the wavelength range.

Acceleration and braking

The PSDs of the contact forces on the front wheelset in straight track portions rolling in acceleration and breaking phases are shown in Figure 25. In the acceleration phases the difference between the experimental and the simulated PSDs is more important on the vertical and longitudinal forces, especially at the shorter wave lengths.

Influence of the number of Karhunen-Loève terms

The numerical simulations with irregularities have been performed with different numbers of considered KLE terms and so different spectral energy relative errors: 10% and 0.1% instead of 1%. The number of needed terms is shown in Table 2 and the evolution of the truncation error according to number of terms is shown in Figure 12. When more terms are taken to divide by 10 the error, no significant changes are observed on the PSDs (Figure 26), but only a little improvement on the vertical force spectrum. In fact, as shown in [START_REF] Panunzio | Wheel-Rail Contact Parameters Sensitivity Analysis with Stochastic Track Geometry[END_REF], the vertical force is mostly sensitive to the superelevation and adding more terms reduces the spectral error on the wavelengths lower than 10 m. The other forces, Y and X are mostly sensitive to the horizontal curvature, which is already well represented on the considered spectral band with a KLE error equal to 1%. When more KLE terms are neglected the spectral energy is lower, especially for the lateral force Y . Indeed, the lateral forces are mostly sensitive to the horizontal curvature irregularities [START_REF] Panunzio | Wheel-Rail Contact Parameters Sensitivity Analysis with Stochastic Track Geometry[END_REF] whose PSD is badly represented for almost all the wavelengths when a KLE truncation equal to 10% is chosen (Figure 8).

Simulation with wheel irregularities

The presence of wheel irregularities can lead to significant contact forces variations. A large wheel flat defect implies large variations on the vertical forces (and their spectrum) because of the impact force induced by a brutal change of the wheel shape [START_REF] Uzzal | Dynamic Analysis of Railway Vehicle-Track Interactions due to Wheel Flat with a Pitch-Plane Vehicle Model[END_REF].

When the wheel untrueness is an inconstant nominal wheel radius (that varies with the wheel rotation angle), periodic changes of the contact forces are produced [START_REF] Newton | An investigation into the dynamics effects on the track of the wheelflats on railway vehicles[END_REF]. This is the situation observed on the experimental spectral densities. Even a small variation on the radius (0.1 mm for instance) induces a variation of the longitudinal creepage (which is normally very small in a straight track at constant running velocity) and so the longitudinal forces. The vertical equilibrium is affected by a very small wheelset roll angle caused by a different rolling radius of the two wheels (approximatively 0.004 • for a 0.1 mm radius difference) that lightly affects the track vertical forces.

Since the wheel nominal diameters of the motor coach (not used for the simulations in this work) are slightly larger (0.46 m) than the passenger car, the peaks appear at shifted wavelengths (Figure 27). Numerical simulations with wheel arbitrary untrueness (nominal radius variable as in Figure 28) on the front wheelset of the leading car are performed and compared with the measurements (Figure 29). The location in frequency of the peaks is well reproduced. Note that the identification of irregularities of the wheels is out of the scope of this paper, so no effort has been devoted to fine tune the amplitude of the peaks.

Conclusion

The steps to lead to a stochastic characterisation of the track geometry irregularities of a railway network have been presented in this work. A large set of experimental measurements is used to model the irregularities as random fields. The model is able to capture the statistical and the spatial variability of the irregularities on the considered track.

Another set of measurements is used to validate the stochastic model of the geometry. The contact forces are measured by another recording train on the same railway network. The stochastic irregularities can be used as inputs of a railway dynamics software to compute the wheel-rail contact forces. The numerical PSDs are very close to the experimental ones when the irregularities are considered in the simulations. The improvement, when comparing the simulations with a theoretical track, is above one order of magnitude and is obtained by considering all the other parameters (such as rail and wheel profiles, friction coefficient and vehicle suspensions stiffness) as deterministic. This means that the track geometry irregularities are the most influential variable parameter.

A further work will consist in investigate on the influence of the track irregularities on the rail rolling contact fatigue, which strongly depends on the wheel-rail contact stresses. As shown in this paper, a deterministic model is not able to capture the real variability of the contact forces (and so the rail surface stresses). The stochastic irregularities model will allow to improve the rail fatigue life predictions.

The method described in the paper only considers the geometry irregularities (i.e. the position of the rails) as stochastic fields modelled from measurements. As shown, when the geometry irregularities are taken into account, a strong improvement of the spectral densities of the contact forces is observed. Other kinds of track irregularities (such as rail corrugation, rail cant variations and joints) could be introduced into the model to further improve the comparison. If these quantities are correlated to the geometry irregularities, then they will be modelled as cross-correlated random fields.

Figure 1 .

 1 Figure 1. Sketch of track geometry irregularities (design position in black, irregular position in green). From top to bottom: gauge, alignment, superelevation, vertical profile.

Figure 2 .

 2 Figure 2. Curvature to versine transfer function (Equation2): (left) module and (right) phase.

Figure 3 .

 3 Figure 3. Example of track geometry measurement with corresponding theoretical values (dashed red line), from left to right: gauge, superelevation, horizontal curvature and vertical curvature.

Figure 4 .

 4 Figure 4. Correlations (at spatial lag equal to zero) between the design horizontal curvature (|C Hd (s)|), the gauge measure (G) and the irregularities fields (gauge G irr , superelevation E irr , horizontal curvature C Hirr and vertical curvature C V irr ).

Figure 5 .

 5 Figure 5. Experimental PDFs of geometry irregularities, from left: gauge, superelevation, horizontal curvature and vertical curvature.

Figure 6 .

 6 Figure 6. Joint PDFs of the track irregularities (at zero lag). Logarithmic color scale.

Figure 7 .

 7 Figure 7. Correlation function of geometry irregularities, from left: gauge, superelevation, horizontal curvature and vertical curvature. Different scales of the lag axis (x).

Figure 8 .

 8 Figure 8. Experimental PSD of geometry irregularities, from left: gauge, superelevation, horizontal curvature and vertical curvature.

Figure 9 .

 9 Figure 9. Correlation function of design horizontal curvature.

Figure 10 .Figure 11 .

 1011 Figure 10. Histogram (normalized) of the length between two welds superposed with the correlation of the horizontal curvature irregularity.

Figure 12 .Figure 13 .

 1213 Figure 12. Karhunen-Loève truncation errors.

Figure 14 .

 14 Figure 14. Percentage error in the PCE identification (ε P CE in Equation12) of the vertical curvature irregularity (100 KLE terms). The dashed red line delimits the area in which the identification is not possible since the number of PCE terms is lower than the number of KLE terms to model. The dotted white line delimits the area in which the PCE error is lower than 1%.

Figure 15 .

 15 Figure[START_REF] Claus | Modeling and Simulation of Railway Bogie Structural Vibrations[END_REF]. PSDs of geometry irregularities, from left: gauge, superelevation, horizontal curvature and vertical curvature. Measurements (red line), generations with 10% (dotted green line with circle markers), 1% (dashed blue line) and 0.1% (dotted violet line with square markers), 1% of the mean spectral density (horizontal dashed black line).

Figure 16 .

 16 Figure[START_REF] Iyengar | Random field modeling of track irregularities[END_REF]. PDFs of geometry irregularities, from left: gauge, superelevation, horizontal curvature and vertical curvature. Measurements (red line), generations (dashed blue line).

Figure 17 .

 17 Figure 17. Example of geometry irregularities, from left: gauge, superelevation, horizontal curvature and vertical curvature. Measurements (red line), generations (blue line).

Figure 18 .

 18 Figure 18. Contact forces decomposition (red) and track reference frame (black).

Figure 19 .

 19 Figure 19. Measurement vehicle schema.

Figure 20 .

 20 Figure 20. Position of the strain gauges used for the contact forces measurement on a wheel.

Figure 21 .

 21 Figure 21. Example of the measured contact forces (front wheelset, front bogies, leading car) with the corresponding horizontal curvature, superelevation and train velocity.

Figure 22 .

 22 Figure 22. Example of contact positions in a curved track

Figure 23 .

 23 Figure[START_REF] Pombo | An alternative method to include track irregularities in railway vehicle dynamic analyses[END_REF]. PSD of contact forces, front wheelset in alignment, constant velocity. Experimental data (red line), numerical simulations with (dashed blue line) and without (dotted black line) irregularities.

Figure 24 .

 24 Figure[START_REF] Luber | Track geometry evaluation method based on vehicle response prediction[END_REF]. PSD of contact forces, front wheelset in curve (1st row outer rail, 2nd row inner rail), constant velocity. Experimental data (red line), numerical simulations with (dashed blue line) and without (dotted black line) irregularities.

Figure 25 .

 25 Figure[START_REF] Li | Assessing track geometry quality based on wavelength spectra and track-vehicle dynamic interaction[END_REF]. PSD of contact forces, front wheelset in alignment, variable velocity (1st row acceleration, 2nd row breaking). Experimental data (red line), numerical simulations with (dashed blue line) and without (dotted black line) irregularities.

Figure 26 .

 26 Figure[START_REF] Kanehara | Measuring rail/wheel contact points of running railway vehicles[END_REF]. PSD of contact forces, front wheelset in alignment, constant velocity. Experimental data (red line), numerical simulations with track geometry irregularities with KLE error equal to 10% (dotted green line with circle markers), 1% (dashed blue line) and 0.1% (dotted violet line with square markers).

Figure 27 .

 27 Figure 27. Experimental PSD of the leading wheelsets of the passenger car (red line) and motor coach (dotted grey line).

Figure 28 .Figure 29 .

 2829 Figure 28. Wheel untrue used for the numerical simulations.

Table 1 .

 1 Number of terms used in the random fields model. KLE truncation error ε KL = 0.01

	December

Table 2 .

 2 Number of terms used in the random fields model with two different KLE truncation error, where d and ng indicate the polynomial chaos degree and the number of germs.

	KLE truncation error		0.1		0.01		0.001
	No. of terms	KLE	PCE (Hermite)	KLE	PCE (Hermite)	KLE	PCE (Hermite)

Appendix B: Continuity of the generated fields

Let l l and l r be respectively the left and the right random limits of the generated field g(s) in Equation 19at the breaking point kl (with k going from 1 to n -1):

By using the expectations in Equation 20 and the definition of the KLE in Equation 5, it follows that, when the number of KLE terms is enough large to give a small truncation error:

As shown in [START_REF] Makhoul | On the eigenvectors of symmetric Toeplitz matrices[END_REF][START_REF] Unser | On the approximation of the discrete Karhunen-Loeve transform for stationary processes[END_REF], the stationarity implies that the KLE modes are alternatively symmetric or anti-symmetric (ϕ i (s) = ±ϕ i (L -s)). Because of this property, it follows that:

By using the definition of the matrix K (Equation 16) and the expectations in Equation 20, it follows that:

From Equations 23 and 24, it follows that the left and right limits are approximate equal:

When M → ∞, the approximate equalities in the Equations 22-25 become exact equalities. Each generated realisation of the random field is thus continuous when M → ∞:

Note that the KLE yields the best such basis in the sense that it minimizes the total mean squared error (Equation 12). The basis is optimal for the decomposition of the covariance function C(|s -t|). The Equation 16 is a decomposition of C(|s -t -L|) into a basis which is optimal for C(|s -t|). For this reason:

This fact also comes from the Cauchy-Schwarz inequality:

December 17, 2016 Vehicle System Dynamics paper In numerical applications, when a large number of KLE terms is chosen, the error made by considering the previous approximate equalities as exact equalities is small. In this work this error is evaluated through the following quantity (r) for the four generated fields:

In Figure 31, the value of this quantity, calculated for different numbers of KLE terms and for each field, is shown. For the numbers of KLE terms chosen for the generation (Table 1), the order of magnitude of r is 10 -2 .

No. of terms [-]