Pierre Aumond

Arnaud Can

Vivien Mallet

Bert De Coensel

Carlos Ribeiro

Dick Botteldooren

Catherine Lavandier

Kriging-based spatial interpolation from measurements for sound level mapping in urban areas

Keywords: Kriging, Spatial interpolation, Sound map, Opportunistic measurements I

come

Introduction

In this work, we investigated the use of a large amount of in situ measurements and the interpolation of these to construct a sound level map.

The Directive 2002/CE/49 contributed to the development and harmonization of noise prediction models (EC 2002). For making urban sound maps, model-based numerical engineering methods are currently widely used and these methods provide a good compromise between accuracy and computation time [START_REF] Kephalopoulos | Advances in the Development of Common Noise Assessment Methods in Europe: The CNOSSOS-EU Framework for Strategic Environmental Noise Mapping[END_REF]. Nevertheless, they have many limitations, and the resulting sound maps neglect the diversity of urban sound environments in terms of both sound sources and sound environment dynamics.

Sound maps based on measurements can help to improve sound mapping [START_REF] Zambon | The LIFE DYNAMAP Project: Towards a Procedure for Dynamic Noise Mapping in Urban Areas[END_REF][START_REF] Asensio | Acoustics in Smart Cities[END_REF][START_REF] Hong | Exploring Spatial Relationships among Soundscape Variables in Urban Areas: A Spatial Statistical Modelling Approach[END_REF][START_REF] Harman | Performance Evaluation of IDW, Kriging and Multiquadric Interpolation Methods in Producing Noise Mapping: A Case Study at the City of Isparta, Turkey[END_REF]. The recent development of small and autonomous acoustic sensors contributes to this movement, and network-based sound monitoring systems are being deployed in an increasing number of cities over the world, using either high-quality or low-cost micropones [START_REF] Mydlarz | The Implementation of Low-Cost Urban Acoustic Monitoring Devices[END_REF][START_REF] Asensio | Acoustics in Smart Cities[END_REF][START_REF] Sevillano | DYNAMAP -Development of Low Cost Sensors Networks for Real Time Noise Mapping[END_REF]. Also, smartphone applications, allowing participatory sensing, are now common, which multiplies the amount of available data to potentially map the sound environment of a city based on measurements [START_REF] Aspuru | CITI-SENSE: Methods and Tools for Empowering Citizens to Observe Acoustic Comfort in Outdoor Public Spaces[END_REF][START_REF] Guillaume | Noise Mapping Based on Participative Measurements[END_REF][START_REF] Maisonneuve | NoiseTube: Measuring and Mapping Noise Pollution with Mobile Phones[END_REF][START_REF] Issarny | Dos and Don'ts in Mobile Phone Sensing Middleware: Learning from a Large-Scale Experiment[END_REF].

Nevertheless, the production of sound maps based on measurements is complicated by metrological issues inherent to the typical microphones used in consumer electronics, the time sparseness of the measurements collected through mobile monitoring applications, and the space sparseness of the measurements collected through fixed sound monitoring networks. Therefore, it is fundamental to know the time and space representativeness of such measurements, as this knowledge is required to be able to propose relevant interpolation methods that can be used to produce sound maps that cover the full temporal and spatial extent of the study area.

The temporal structure of urban sound levels (highly correlated day or week patterns, seasonal trends) can be exploited to restrict the number of sampled days, or to rely only on measurements performed at selected periods of the day, in order to estimate Lden values or Daily Average Noise Patterns [START_REF] Hong | Exploring Spatial Relationships among Soundscape Variables in Urban Areas: A Spatial Statistical Modelling Approach[END_REF][START_REF] Geraghty | Investigating the Temporal Variability of Noise in an Urban Environment[END_REF]; J. M. [START_REF] Barrigón Morillas | Uncertainty Evaluation of Continuous Noise Sampling[END_REF][START_REF] Zuo | Temporal and Spatial Variability of Traffic-Related Noise in the City of Toronto, Canada[END_REF]. Previous studies also revealed that a 10 or 15-minute measurement is representative of a 1h-period in an urban context, as the majority of the 10 or 15 minute-measurements are in the same sound level range during homogeneous periods (Brocolini et al. 2013;[START_REF] Prieto Gajardo | Stabilisation Patterns of Hourly Urban Sound Levels[END_REF]. Other studies that rely on short-term recording methods proved the relevance of 15-minute sampling periods [START_REF] Morillas | A Categorization Method Applied to the Study of Urban Road Traffic Noise[END_REF][START_REF] Can | Sampling Approaches to Predict Urban Street Noise Levels Using Fixed and Temporary Microphones[END_REF]. Even shorter measurement periods can be found in the literature, especially in those cases where the opportunistic measurement context offered by smartphone applications is used. In this case, the short measurement duration is compensated by the large number of measurements [START_REF] Guillaume | Noise Mapping Based on Participative Measurements[END_REF], shifting the focus from the duration of each measurement episode to the number of sampling episodes, as recommended in [START_REF] Mateus | Assessing the Influence of the Sampling Strategy on the Uncertainty of Environmental Noise Measurements through the Bootstrap Method[END_REF]. [START_REF] Zuo | Temporal and Spatial Variability of Traffic-Related Noise in the City of Toronto, Canada[END_REF] showed that the sound level variability in urban environments can be explained for a large part by the spatial characteristics of the environment. Also, the space representativeness and the spatial interpolation of the measurements is an important issue when computing sound maps based on measurements. Maps interpolated from the data obtained through fixed sound measurement stations have recently been produced [START_REF] Liu | Spatiotemporal Variability of Soundscapes in a Multiple Functional Urban Area[END_REF][START_REF] Harman | Performance Evaluation of IDW, Kriging and Multiquadric Interpolation Methods in Producing Noise Mapping: A Case Study at the City of Isparta, Turkey[END_REF][START_REF] Garcia | Spatial Statistical Analysis of Urban Noise Data from a WASN Gathered by an IoT System: Application to a Small City[END_REF][START_REF] Huang | Optimizing Stations Location for Urban Noise Continuous Intelligent Monitoring[END_REF], and form a useful tool to estimate the noisiness of a neighborhood or to give a global overview of the city sound levels. However, the large distance between measurement stations often does not allow one to map sound levels in each street, which is offered by maps calculated using model-based numerical methods. A study by [START_REF] Can | Measurement Network for Urban Noise Assessment: Comparison of Mobile Measurements and Spatial Interpolation Approaches[END_REF] that involved mobile measurements performed using sound level meters attached to backpacks permitted one to compare an interpolated map with a reference map, but only for a very small area (four streets). More studies are therefore needed to investigate the density of measurements that is required to have an acceptable accuracy at the street resolution.

Another parameter to take into account is the method of interpolation. Several methods have been tested for urban sound level interpolation: Inverse Distance Weighting (IDW) methods (A. [START_REF] Can | Measurement Network for Urban Noise Assessment: Comparison of Mobile Measurements and Spatial Interpolation Approaches[END_REF][START_REF] Harman | Performance Evaluation of IDW, Kriging and Multiquadric Interpolation Methods in Producing Noise Mapping: A Case Study at the City of Isparta, Turkey[END_REF][START_REF] Garcia | Spatial Statistical Analysis of Urban Noise Data from a WASN Gathered by an IoT System: Application to a Small City[END_REF][START_REF] Hong | Exploring Spatial Relationships among Soundscape Variables in Urban Areas: A Spatial Statistical Modelling Approach[END_REF], Kriging methods [START_REF] Harman | Performance Evaluation of IDW, Kriging and Multiquadric Interpolation Methods in Producing Noise Mapping: A Case Study at the City of Isparta, Turkey[END_REF][START_REF] Can | Measurement Network for Urban Noise Assessment: Comparison of Mobile Measurements and Spatial Interpolation Approaches[END_REF][START_REF] Garcia | Spatial Statistical Analysis of Urban Noise Data from a WASN Gathered by an IoT System: Application to a Small City[END_REF]) and multiquadratic interpolation [START_REF] Harman | Performance Evaluation of IDW, Kriging and Multiquadric Interpolation Methods in Producing Noise Mapping: A Case Study at the City of Isparta, Turkey[END_REF]. Alternative interpolation methods that involve a modified definition of the distance between measurement locations, in order to account for the city geometry or the road network, have recently been proposed (A. [START_REF] Can | Measurement Network for Urban Noise Assessment: Comparison of Mobile Measurements and Spatial Interpolation Approaches[END_REF][START_REF] López-Quílez | Geostatistical Computing of Acoustic Maps in the Presence of Barriers[END_REF][START_REF] Hachem | Monitoring Noise Pollution Using the Urban Civics Middleware[END_REF]. Nevertheless, these interpolation methods have only been tested on small measurement samples, and larger studies are needed to validate the conclusions of these works. It is worth mentioning that the density of measurements appears to be more important than the method of interpolation [START_REF] Harman | Performance Evaluation of IDW, Kriging and Multiquadric Interpolation Methods in Producing Noise Mapping: A Case Study at the City of Isparta, Turkey[END_REF].

Recently, a number of data fusion techniques have been proposed to correct model-based numerically computed sound maps with measurements [START_REF] Hachem | Monitoring Noise Pollution Using the Urban Civics Middleware[END_REF][START_REF] Wei | Dynamic Noise Mapping: A Map-Based Interpolation between Noise Measurements with High Temporal Resolution[END_REF][START_REF] Zambon | The LIFE DYNAMAP Project: Towards a Procedure for Dynamic Noise Mapping in Urban Areas[END_REF][START_REF] Ventura | Estimation of Urban Noise with the Assimilation of Observations Crowdsensed by the Mobile Application Ambiciti[END_REF]). These techniques are promising, but require a pre-calculated sound map, which can be expensive, and give priority to the potential indicators available in this sound map. Most of the existing sound maps involve the energy-equivalent sound level, whereas in situ measurements allow calculating a wide range of acoustical indicators, which may include information about the temporal dynamics of the sound environment; for example, percentile or emergence indicators can also be interesting to interpolate. The spatial interpolation methods can even be based on perceptual assessments [START_REF] Aletta | Soundscape Approach Integrating Noise Mapping Techniques: A Case Study in Brighton, UK[END_REF].

In this study, a large measurement campaign has been conducted in the XIII th district of Paris, with the goal to test different spatial interpolation strategies. Mobile measurements have been performed with sound measurement stations attached to backpacks that were carried by researchers when walking in every street of the district between 1 and 15 times. The measurements are aggregated over a grid of locations in the study area, and are used to compute a reference map of the district. An analysis of the sensitivity of the sound level values with respect to the radius of the integration and the number of measurements is done via a bootstrap method. From the reference sound map, four Kriging methods for interpolation between a set of measurement locations are tested, based on a combination of two strategies: (i) Ordinary Kriging and universal Kriging which consists of adding a linear trend, defined from the distance between an observation location of the domain and its closest categorized road, and (ii) a variation of the distance definition between observation locations, which can be Euclidian or computed from the road network to take into account the influence of the city geometry. By progressively decreasing the number of observation locations, the impact of the density of observation locations and the performance of different spatial interpolation methods is investigated.

II. Method

A. Study area

Figure 1 presents the study area, which corresponds to the XIII th District of Paris. This district includes a large variety of urban sound environments: large avenues with high traffic density, lively streets with bars and restaurants, schools, small and large parks, quiet streets. The size of the study area is approximately 2.8 km² with a maximum extent of 2 km west to east and a maximum extent of 1.7 km north to south.

B. Measurement set-up

All measurements were carried out using dedicated mobile sound measurement stations developed by ASAsense [START_REF] Coensel | Dynamic Noise Mapping Based on Fixed and Mobile Sound Measurements[END_REF]. These stations record the instantaneous 1/3-octave band spectrum with a 125-ms temporal resolution as well as the instantaneous GPS position with a 1-s temporal resolution. Both sound and GPS data are synchronized, such that the spatio-temporal evolution of the sound spectrum during each measurement session can be reconstructed afterwards. In order to fully capture the characteristics of the sound environment, a large set of indicators is calculated on the basis of this 1/3-octave band spectrum data, including the A-weighted energy-equivalent sound level at each second, L eq,1s , which is used throughout this study.

C. Mobile measurements

The sound measurement stations were mounted in backpacks with power provided by a battery pack, and, subsequently, mobile measurements were carried out between October 22 th 2014 and May 26 th 2015. Five operators participated in the measurements. In order to minimize the variation between measurement sessions and to be able to calculate sound levels that are representative for homogeneous sound environments, measurements were only carried out on weekdays, either between 10 a.m. and 12 a.m., or between 2 p.m. and 4 p.m. As shown in (Prieto Gajardo and Barrigón Morillas 2015;[START_REF] Zuo | Temporal and Spatial Variability of Traffic-Related Noise in the City of Toronto, Canada[END_REF]Brocolini et al. 2013), these periods, which exclude rush hour traffic and lunch times, provide a similar sound environment. Depending on the variability of the sound environment, the number of walks in each street was varied. As shown in (Prieto Gajardo and Barrigón Morillas 2015), the sound environment of a calm street is more sensitive to single sound events than that of a large boulevard, thus a higher number of measurements is required to calculate a sound level value that is representative of the sound environment of a calm street, as compared to a large boulevard. After each day of measurements, the variance of the sound level was computed for each street, providing feedback on those streets that would benefit from more measurements to get a stable estimate. The number of passages per street ranged from 1 to 15 times, with an average of 4 and a standard deviation of 2.7 passages. Figure 1 shows the number of walks (only validated measurements were kept to plot this figure).

D. Map matching

A GPS track is associated with each measurement session. However, the accuracy of the GPS data depends on many factors such as the quality of the GPS receiver, the characteristics of the surroundings such as the presence of high buildings, or the weather. In this study, the median standard deviation associated with the GPS locations was about 10 meters. Although rarely considered in studies dealing with geo-referenced mobile measurements, this can be problematic for the present analysis, because some measurements can be associated with an erroneous street.

Therefore, it was necessary to preprocess the GPS data, mapping each measurement to a location on the road, a problem commonly known as map-matching. For this study, a point-to-point method was developed, on the basis of the following conditions: (i) all measurements were performed on sidewalks or on roads, and the sound levels measured on both sides of each street are considered to be equivalent and are snapped to the middle of the street; (ii) all GPS locations are snapped to the center of the closest street under the conditions that the operator walked with a maximum speed of 5 km/h and that the map-matched point conserves the same direction of displacement (with a direction tolerance of 60°) as the original point; and (iii) the map-matched points are located at a maximum distance from the original GPS points equal to twice the standard deviation of the GPS tracker.

E. From mobile measurements to observation locations

In a previous study [START_REF] Aumond | Modeling Soundscape Pleasantness Using Perceptual Assessments and Acoustic Measurements Along Paths in Urban Context[END_REF], it was shown that the median sound level L 50 is well correlated with the perceived loudness of an urban sound environment. The median sound level presents the advantage that it is less sensitive to peaks in the measurement than the energy-equivalent level. Exceptionally, peaks can be generated by the operators themselves, or can occur when extremely noisy vehicles pass by in the vicinity of the operators. In addition, the median sound level L 50 does not include A-weighting, which is known to reduce too much the influence of low frequencies [63-500 Hz] at sound levels encountered in urban environments, and thus the influence of road traffic sound on overall perceived loudness. The first step in carrying out the spatial interpolation between measurement locations, was to aggregate, at each location, all mobile measurements that are within a radius r. This aggregation step is performed using the median sound level of all the 1-s values L eq,1s considered as independent observations. The associated value is assumed to be representative of the L 50 sound level at the observation location on weekdays and during the measurement periods [10-12h; 14-16h]. Every aggregated observation location is situated on the road network, because mobile measurements where only taken on sidewalks and pedestrian walkways in the public space. The road network that was used for mapmatching was based on OpenStreetMap ("OpenStreetMap" n.d.).

F.

Variogram and Kriging

Kriging method

Ordinary Kriging method is a well-known interpolation method that has been used in various applications, especially in environmental applications. It also bears similarities with classical data assimilation methods that have helped environmental forecasting, including at urban scale for air pollution [START_REF] Tilloy | BLUE-Based NO2 Data Assimilation at Urban Scale[END_REF]) and noise pollution [START_REF] Ventura | Estimation of Urban Noise with the Assimilation of Observations Crowdsensed by the Mobile Application Ambiciti[END_REF][START_REF] Harman | Performance Evaluation of IDW, Kriging and Multiquadric Interpolation Methods in Producing Noise Mapping: A Case Study at the City of Isparta, Turkey[END_REF]. The approach is likely to succeed when some meaningful function can fit the empirical variogram, which is the case in this study. Nonetheless, at urban scale, the sound levels exhibit a special distribution due to the city geometry, and this is a difference with typical applications of Kriging [START_REF] López-Quílez | Geostatistical Computing of Acoustic Maps in the Presence of Barriers[END_REF]. This is the reason why, in the construction of the variogram, we made use of the distance along the road network instead of a raw Euclidean distance (see Section II.F.3). Also, since it is known that sound levels can be approximated from external data (Juan Miguel Barrigón [START_REF] Morillas | Analysis of the Prediction Capacity of a Categorization Method for Urban Noise Assessment[END_REF]) and can provide some sort of prior for the spatial distribution of the sound levels, we also tested universal Kriging, which carries out the interpolation of the measurements on top of a prior sound level map (see Section II.F.4).

Implementations and parameters

The variogram and Kriging algorithms presented in this study are applied using the functions variog (computation of the variogram), variofit (best fit of the variogram) and krige.conv (Kriging function) of the packages GeoR ("GeoR: Analysis of Geostatistical Data Version 1.7-5.2 from CRAN" n.d.) and GeoRcb (López-Quílez and Muñoz 2009). The empirical variogram is computed over a distance of 1000 meters. The classical estimator is chosen to compute the empirical variogram as defined in [START_REF] Cressie | Statistics for Spatial Data[END_REF]. The Matérn covariance model, as defined in [START_REF] Diggle | Directive 2002/49/EC of the European Parliament and the Council of 25 June 2002 Relating to the Assessment and Management of Environmental Noise[END_REF], is used to compute the best fit of the experimental variogram with a fixed value for the shape parameter κ=0.5 and φ the value of the range parameter to estimate.

Euclidian vs cost-based distance

The package GeoR has been used to interpolate using the Euclidian distance (EUCL). The package GeoRcb permits the use of alternative definitions of distance. In this study, the distance between two observation locations has also been defined along the road network as presented in Figure 2 The error correlations in the traffic flow are assumed to be better modeled as a function of the distance along the road than the Euclidean distance. As a consequence, the use of the distance along the road network presumably allows one to better model errors that come from the traffic. The distance between all the possible observation locations of the domain has been calculated with the Johnson's algorithm of the "distances" function of the package "igraph" of the R software described in [START_REF] West | Introduction to Graph Theory -Second Edition[END_REF].

Ordinary and universal Kriging

Two Kriging methods are compared: Ordinary Kriging (OK) and Universal Kriging (UK). Universal Kriging is a variant of the ordinary Kriging operation that includes a linear trend. Barrigon et al. showed that urban sound is strongly stratified [START_REF] Morillas | A Categorization Method Applied to the Study of Urban Road Traffic Noise[END_REF]. Based on a similar statement, the linear trend T in this study is defined as a linear regression based on 4 explanatory variables which are the distances D i between the observation locations of the study area and its closest roads belonging to 4 categories i. Road categories have been defined based on the OpenStreetMap attributes as shown in Table 1. Figure 3 presents a map of the study area, in which the roads are shown in different colors according to their category. The first three categories correspond to streets with vehicular traffic, the last one to pedestrian streets.

Open Street Map attributes

Categories

G. Performance metrics

The performance of the interpolation methods is assessed with two indicators: the Root Mean Square Error (RMSE) and the Pearson correlation coefficient (r) between the interpolated and the reference map. The geospatial interpolations are proposed following four methods: (a) ordinary Kriging (OK + EUCL), (b) ordinary Kriging using the distance along the road network (OK + ROAD), (c) universal Kriging (UK + EUCL), (d) universal Kriging using the distance along the road network (UK + ROAD).

III. Results

A. Reference sound map

A reference sound map based on all mobile measurements is first computed. For each observation location, the mobile measurements that are within a specific integration radius are mapped to this location. A sensitivity analysis of the spatial representativeness and the expected accuracy of the sound level at the aggregation radius and the aggregated number of 1-s samples is then performed.

Six values for the aggregation radius have been tested: 2.5 m, 5 m, 10 m, 15 m, 25 m and 50 m. For each radius value, a set of observation locations has been selected, uniformly distributed. For this statistical analysis, the observation locations are spaced apart by at least two times the radius value, to avoid redundancy. In Figure 4, three of those subsets of observation locations are shown. Figure 5(a) shows the distribution of the aggregated sound levels L 50 over the study area for the smallest and largest aggregation radii (2.5m and 50m) as presented in Section II.E. Figure 5 (b) shows the P 90 -P 10 indicator (where P 10 and P 90 are the percentiles 10 and 90 of the sound level distributions), which reflects the width of these distributions for all the studied aggregation radii. As expected, the range of the aggregated sound levels L 50 over the study area decreases as the value of the radius increases. More prosaically, the reference sound map will appear to be blurred. A small radius will give rise to a more detailed reference sound map, but comes with a decreased number of measurements at each location that may therefore no longer be representative. The results show that the influence of the aggregation radius over the sound level distribution is relatively small, with P 90 -P 10 decreasing from 13.5 dB to 10.5 dB. of the distribution of the sound levels (P90-P10) for 6 different aggregation radii A bootstrapping method (bootstrp function, statistical toolbox, Matlab) is proposed to analyze the sensitivity of the L 50 value to the aggregation radius and the number of measurements. This method relies on random sampling, with replacement, of the 1-s measurements for each location within the study area [START_REF] Efron | Bootstrap Methods: Another Look at the Jackknife[END_REF]. Multiple replications of the method permits the computation of the variance associated with the average L 50 value at each observation location due to the sample characteristics. Figure 6 presents (a) the standard deviation of 1000 bootstrap replications of the calculated L 50 varying with the number of 1-s samples and (b) the relationship between the aggregation radius and the proportion of retained locations, considering a minimum number of 1-second samples of 180. The correlation between these two parameters is statistically significant (r=0.53, p<.05), but a large aggregation radius does not always imply a large number of measurements (e.g., along the borders of the study area) and vice versa (e.g., if the operator measured a few minutes at one specific location). The variation specific to each location is found to be considerably smaller than the global variation between the locations over the whole study area, which is about 5 dB (see Figure 5). If we consider a standard deviation smaller than 1 dB as acceptable, 90% of the calculated standard deviations are below this threshold when a sample is composed by a minimum of 180 1-s measurements. The chosen threshold also guarantees to have enough data to carry out the spatial analysis. An aggregation radius of 25 meters is chosen as the threshold value because it permits one to retain more than 60% of the computed locations. This aggregation radius also has the advantage that it corresponds to the longitudinal spatial representativeness of sound level measurements found in literature, perceptually [START_REF] Brocolini | Discrimination of Urban Soundscape through Kohonen Map[END_REF] or physically (A. [START_REF] Can | Measurement Network for Urban Noise Assessment: Comparison of Mobile Measurements and Spatial Interpolation Approaches[END_REF], thus suggesting that the resulting sound level map will be consistent with observed spatio-temporal variations.

On the basis of the above results, it can be concluded that, at least for the urban study area considered, three-minute measurements provide sufficient confidence in the aggregated measurement value, even if its representativeness is not guaranteed. As the purpose of the present study was to perform a comparison of interpolation techniques, the calculated sound maps only need to respect the spatial variation of the sound level, but do not necessarily need to be representative of a homogeneous time period.

Figure 7 presents the resulting reference map of the median sound level (L 50 , dB). For the calculation of this map, a distance of 10 m between each observation location was used, resulting in 4360 locations.

B. Spatial dependence of the data

The spatial dependence of the data is highlighted through the calculation of variograms, which express the semivariance between L 50 values for a couple of locations according to their distance. On Figure 8 The parameters of the best fitted covariance models are presented in Table 1. The practical ranges of the variograms, defined as the value for which the correlation function decays to 5% of its value at 0, are (a) 366 m for OK+EUCL (b) 285 m for UK+EUCL (c) 691 m for OK+ ROAD (d) 481 m for UK+ROAD. For OK+EUCL, hardly any information is given by an observation to an estimated value located in a radius superior to 366 m. As the nugget variance τ² is null, the variogram asymptote, or sill value, corresponds to the signal variance σ². This value at 1000 meters is 22.4 dB² and 24.2 dB² for OK + EUCL and OK + ROAD methods, and 11.2 dB² and 11.9 dB² for UK + EUCL and UK + ROAD. Thus, adding the trend permits one to considerably reduce the signal variance and illustrates the strong correlation between the urban sound levels and the proximity to different types of roads. Also, the practical range of the variograms increases when the alternative definition of distance from the road network is used (from 366 to 691 m and from 285 to 481 m).

Spatial interpolation and performance analysis

A subset of the observation locations is selected from the reference map and is interpolated over the whole study area using the four tested strategies. As described in Section III.A, the median sound level values at the observation locations are computed from at least 180 L eq,1s measurements included in a 25 m radius around each observation location. The subset of observation locations for evaluation is randomly selected and 1000 replications are performed for each subset configuration. The replications give information about the variability due to a chosen set of observation nodes. Figure 9 presents the relationships between two indicators of performance (RMSE and the Pearson correlation coefficient), and the density of nodes (number of nodes per sq. km) for the four methods of interpolation. Universal Kriging (UK) considerably increases the quality of the results on both indicators compared to ordinary Kriging (OK). Taking into account the distance along the road network (ROAD) only leads to better results for the ordinary Kriging cases.

From a practical point of view, in case only fixed measurement stations are used, 15 observation locations per sq. km is already a large number. In this case, Figure 9 shows that the correlation between interpolated and reference sound levels is between 0.5 and 0.8, and that the RMSE value is between 2.5 and 3.5 dB. Figure 10 shows an interpolated sound map, using the UK + EUCL and the OK + EUCL methods, based on one of the random sets of 42 observation locations (15 locations per sq. km). The associated prediction standard error maps are also given, as the Kriging methods give access to this information [START_REF] Diggle | Directive 2002/49/EC of the European Parliament and the Council of 25 June 2002 Relating to the Assessment and Management of Environmental Noise[END_REF]. Figure 10 shows a good correspondence with the reference sound map (see Figure 7). Nevertheless, the dispersion in Figure 9 shows that the observation locations have an important influence over the performance of the algorithm. Figure 11 shows the average errors for the 1000 replications (42 observation locations), against the reference map presented in Figure 7. Figure 11 (a-d) shows that all the interpolation strategies underestimate the sound level for large boulevards and overestimate the sound level for quiet places, but this is a common limitation of interpolation methods. This is less pronounced for universal Kriging (Figure 11 Figure 9 shows that the standard deviation of the performance indicators can be important, especially in case of a small density of measurement locations. In order to identify the principal factors which influence the performance of the Kriging algorithms, a short statistical analysis has been performed over geospatial indicators of the spatial distribution of 42 observations. Table 2 presents some of the parameters which have been calculated for each generated set of observations. The Variance Mean Ratio (VMR), also called index of dispersion, is an indicator of a good dispersion of the observations over the domain. For this, the domain is divided in 6x6 cells. In each cell, the number of observations is calculated. The variance mean ratio corresponds to the ratio between (i) the variance of the number of observations between each cell, and (ii) the mean number of observations for a cell. If VMR is equal to 1, the observations are randomly dispersed; if VMR<1 the observations are more dispersed than random (e.g., regular distribution); if VMR > 1 the observations are more clustered than a random dispersion. MEAN_ROAD, STD_ROAD, KURT_ROAD, SKEW_ROAD are indicators of the distribution of the observations over the road categories (integers from 1 to 4, see Figure 2). For example, if MEAN_ROAD is equal to 1 and STD_ROAD equal to 0, it means that all the observations are located over roads of type 1 (avenue/boulevard).

Finally, the CENTER indicator represents the distance between the center of gravity of the observation locations and the center of gravity of the whole road network within the case study area. If CENTER is small, the observations are well centered over the domain. Table 3 presents the Pearson correlation coefficient between the spatial distribution indicators and the performance indicators. It shows that the STD_L50 is the only geospatial indicator that is well correlated with the performance indicators, particularly for ordinary Kriging. As it can be expected, this suggests that the measurement stations should be placed in locations that represent a large distribution of sound levels. For universal Kriging, this factor is less important because of the correction brought by the trends. The same analysis was carried out with other performance metrics than the Pearson correlation coefficient, but this did not bring more information.

IV. Discussion

As the study area was the XIII th district of Paris, the variogram and Kriging parameters were only adjusted for this district. Even if some of these parameters are comparable with those from previous studies (A. [START_REF] Can | Measurement Network for Urban Noise Assessment: Comparison of Mobile Measurements and Spatial Interpolation Approaches[END_REF], other replications in other urban contexts will be necessary to extend the conclusions of this work.

On the one hand, the alternative definition of distance along the road network slightly increases the performance of the algorithms, but only for ordinary Kriging methods. On the other hand, the proposed trend based on the distance to the closest roads by category, as defined in Section II.F.4, strongly improves the results for all Kriging methods. These results confirm the observations done in [START_REF] Morillas | A Categorization Method Applied to the Study of Urban Road Traffic Noise[END_REF][START_REF] Gozalo | Urban Noise Functional Stratification for Estimating Average Annual Sound Level[END_REF][START_REF] Morillas | Analysis of the Prediction Capacity of a Categorization Method for Urban Noise Assessment[END_REF] which state that using a street categorization method that accounts for street use is particularly appropriate to study the spatial variability of urban sound levels. Thus, to continue to improve the performance of the interpolation methods, efforts should be focused on the trend definition. Nevertheless, a complex trend definition is comparable to the computation of a model-based sound map. Moreover, it might increase the computational cost of the interpolation method considerably, and can possibly introduce new error sources.

The errors associated with the observations are not taken into account in our Kriging approach; nevertheless, a follow-up study should investigate the integration of noisy observations with a proper uncertainty associated with each measurement location. Maybe data assimilation methods used in geosciences, e.g., computing the so-called best linear unbiased method [START_REF] Bouttier | Data Assimilation Concepts and Methods[END_REF][START_REF] Tilloy | BLUE-Based NO2 Data Assimilation at Urban Scale[END_REF], could be useful to achieve this task.

All the presented results are valid for one homogeneous time period. Time interpolation could be added. This interpolation could be done in pre-or post-processing, relying on previous works dealing with temporal interpolation of urban sound levels (Prieto [START_REF] Gajardo | Can Weekly Noise Levels of Urban Road Traffic, as Predominant Noise Source, Estimate Annual Ones?[END_REF] or directly with a spatio-temporal Kriging.

Finally, this study shows the influence of the spread of observation locations over the study area to correctly interpolate sound levels. As it can be expected, the results suggest that fixed measurement stations should be placed to obtain a large distribution of sound levels and so, a large variety of sound environments. Nevertheless, even if it can be estimated, this information is not fully available prior to the installation of the measurement stations. Other indicators, such as the distribution of the measurement stations over the various road categories, or the location of the center of gravity of the measurement stations, which are available a priori, do not show relevant relationship with performance indicators.

V. Conclusion

By means of a progressive degradation of a reference map of 2.8 km² interpolated from geo-referenced mobile measurements, spatial interpolation methods were compared. The impact of the density of observation points and the performance of four spatial interpolation methods were presented. The four interpolation methods were constructed by combining two algorithms: (i) the Kriging method, either ordinary Kriging or universal Kriging (which consists in adding a linear trend, defined from the distance between each location and its closest road in each category) and (ii) the definition of the distance between locations, either Euclidian or computed from the road network.

The main conclusions are:

 A minimum of 180 1-s measurements are needed to obtain an acceptable level of confidence (1dB) in the L 50 value calculated at each location within the study area, with an aggregation radius of 25 m.

 The practical ranges of the variograms computed for the four Kriging methods are between 250 and 700 meters.

 Using the distance along the road network in the Kriging method considerably increases the performance in case of ordinary Kriging, but not in case of universal Kriging.

 Universal Kriging, which consists of adding a local trend in the ordinary Kriging formulation, is a promising method. Nevertheless, it introduces an additional calculation of the trend that has a pre-processing cost and can in itself be a source of error.

 Approximately 50 observation locations per sq. km are needed in order to get a correlation coefficient superior to 0.8 and a RMSE value inferior to 2.5 dB between the reference and the interpolated map.

In view of the large density of observation locations needed to obtain a sound map with a high accuracy and the strong improvement brought by the trend in the Kriging formulation, further studies should probably focus on fusion or assimilation techniques combining measurements and numerical simulation results.

Figure 1 (

 1 Figure 1 (color online) Number of passages at each location (only validated measurements).

Figure 2 (

 2 Figure 2(color online) Schema of the road network (shown in black), with the Euclidian Distance (shown in red) as well as the Road Network Distance (shown in green).

  Equation 1 presents the equation of the trend: T ~ a.D 1 +b.D 2 +c.D 3 +d.D 4 +e in which a, b, c, d and e are adjusted constants.

  color online) Outline of the study area (black dashed line) and, in color, the 4 roads categories.

  of three selected sets of observation locations (dots) for radii: (a) 2.5m, (b) 15 m and (c) 50m. Observation locations are located on the study area road network (solid lines).

Figure 5 (

 5 Figure 5 (color online) (a) Distribution of the aggregated sound levels L50 over the study area for two values of the aggregation radius. (b) Width

Figure 6 (

 6 Figure 6 (color online) The upper plot shows the relationship between the standard deviation of 1000 replications (subsets) of the underlying reference data set and the number of 1-s samples of the subsets. The lower plot shows the relationship between the aggregation radius and the proportion of retained measurement locations, considering a minimum number of 1-second samples of 180. Each dot represents a replication, the solid line denotes the mean value, the dashed line denotes the percentile-90).

Figure 7 (

 7 Figure 7 (color online) Reference sound map, integration radius 25 meters, integration time > 180 sec. (n= 4360)

  , four fitted variograms derived from the reference sound map are presented: (a) an ordinary variogram with the Euclidean distance (OK+EUCL), (b) an ordinary variogram using the distance along the road network (OK+ROAD) (c) an universal variogram which accounts for the trend with the Euclidean distance (UK+EUCL), (d) an universal variogram which accounts for the trend and also uses the distance along the road network (UK+ROAD).

Figure 8

 8 Figure 8 Empirical variograms (dots) and best fitted parametrical models (red line) computed using the ordinary Kriging (OK) and universal Kriging (UK) methods. The printed distance is computed using the Euclidian distance (EUCL) or the distance along the road network (ROAD) (distance in meters, and semivariance in dB²).

  Figure 9 (color online) Relationships between the two indicators of quality (RMSE and the Pearson correlation coefficient), and the density of nodes (number of nodes per sq. km) for the four methods of interpolation. The area corresponds to the standard deviation associated with the 1000 replications. (The red and blue lines are nearly superposed).

  Figure 10 (color online) Example of an interpolated sound map from 42 observations: Estimated L50 values using the (a) OK+EUCL and (b) UK+EUCL Kriging methods, and the associated standard errors maps (c) OK+EUCL and (d) UK+EUCL

  Figure 11 (color online) Average errors for 1000 replications (42 observation locations) against the reference map REF for the four interpolation methods (a) OK+EUCL, (b) UK+EUCL, (c) OK+ROAD and (d) UK+ROAD.

Table 1 -

 1 Parameters of the Kriging methods

			Nugget	Sill	Range	Practical Range
	Kriging method	Covariance model	(dB²) variance	(dB²)	(m) parameter	(m)
			τ²	σ²	φ
	OK+EUCL	Matérn with fixed κ = 0.5	0	22.4	122.2	366
	UK+EUCL	Matérn with fixed κ = 0.5	0	11.2	95.1	285
	OK+ROAD	Matérn with fixed κ = 0.5	0	24.2	230.7	691
	UK+ROAD	Matérn with fixed κ = 0.5	0	11.9	160.6	481

Table 2 -

 2 Geospatial indicators of the observation locations.

	Index	Name

Table 3 -

 3 Pearson correlation coefficient between the geospatial and performance indicators for two Kriging

		methods (* p<0.05).
		r	r	RMSE	RMSE
		(OK+EUCL)	(UK+EUCL)	(OK+EUCL)	(UK+EUCL)
	VMR	0.06*		-0.07*
	STD_L50	0.48*	0.27*	-0.55*	-0.21*
	MEAN_ROAD	0.07*			-0.04*
	STD_ROAD	0.19*	0.07*	-0.19*	-0.10*
	KURT_ROAD	-0.10*	-0.06*	0.07*	0.09*
	SKEW_ROAD	-0.16*	-0.07*	0.14*	0.11*
	CENTER	-0.05*		0.06*

Acknowledgments

This work was carried out in the framework of the GRAFIC project, supported by the French Environment and Energy Management Agency (ADEME) under contract No. 1317C0028. The authors also would like to thank the BruitParif team which completed most of the measurements.