Comparison of Stochastic and Deterministic Methods for Mapping Environmental Noise from Opencast Quarries

Corresponding Author: Dario Lippiello Department of Engineering, University of Roma Tre, Rome, Italy Email: dario.lippiello@uniroma3.it Abstract: An accurate mapping of the acoustic spatial variability generated by quarrying plants constitutes an important step in the assessment of the environmental compatibility of the extraction site with the surrounding area. This paper refers to the results of a study comparing performances obtained by stochastic and deterministic methods to determine the acoustic climate of the areas surrounding opencast quarries. To this aim noise levels are monitored in the areas surrounding a limestone quarry site in central Italy. In the first phase an analytical prediction model, based on ISO 9613 and ISO 3746, is developed using the sampled sound power level of each source as input data. In the second phase, in contrast, the stochastic method of Ordinary Kriging is used to plot an acoustic map. This is drawn up by interpolating the scattered sound pressure level data from sampled points which are identified as nodes of a rectangular grid covering the chosen area. The analytical prediction model output, which consists in the sound pressure level at various target points, is then compared with the homologous map obtained from the Ordinary Kriging method. Finally the resulting maps are evaluated with respect to statistical measures and cross validation procedures in order to determine the most suitable method.


Introduction
Quarrying activities and extractive sites, involving both labour and material resources, represent a productive sector which is of basic importance to the Italian economy. These activities have a considerable impact on many aspects of the environment (i.e., Ozcan et al., 2012). One such aspect is that of change in the acoustic climate of the surrounding areas, constituting a significant disturbance for anybody living in the vicinity of the mine or quarry. The effects, however, are not only determined by the plant itself and there is a considerable impact generated by all the activities associated with the plant such as the increase in haulage traffic to and from the site, or the use of explosives for the extraction of raw materials (Giraudi et al., 2009), airborne dust emissions due to quarrying activities (Bluvshtein et al., 2011), rather than visual impact on the surrounding areas (Alfaro Degan et al., 2014).
Historically, the need for acoustic maps, which urgently arose for the interested parties or stakeholders (the public authorities, business owners, private citizens etc.) has been tackled using a number of approaches (Yepes et al., 2009).
The first and most common approach regards the application of the source models (Alfaro Degan et al., 2005) based initially on the characterization of the sound power level generated by the identified sources of the plant (Roy and Adhikari, 2007) and subsequently able to represent the evolution in space of the signal, by means of the application of a physical model. On the basis of ISO 9613 standard, many studies have focused on how best to determine the physical parameters to compute in the sound propagation model (Cinar and Sensogut, 2009;Dowd and Li, 2000), that is on those parameters linked to the geomorphological characteristics of the sites, from those regarding the plant operation (Pathak et al., 1999), to the attenuation effects and changes in acoustic signal due to differences in the types of materials extracted (Neto et al., 2012).
If on the one hand this modelling approach has improved the quality of the results, on the other hand it has brought about an increase in the number of parameters to determine and the measurements to be taken, augmenting the computational complexity of the model as well as increasing the costs and time required to carry out research.
In order to obviate this tendency, new alternative approaches to the physical model are being studied (Nanda et al., 2009). The aim of this article is to carry out and test, with regard to the case in question, a geostatistical answer to the problem of acoustic mapping in the areas surrounding the extractive sites. The sound pressure level is thus considered as a regionalized variable and its spatial distribution is generated by a stochastic process (Baume et al., 2008). The technique utilised here, which is widely used in environmental science, is that of ordinary Kriging. The geostatistical results may then be compared with the measurements taken in the field.

Materials and Methods
This research was carried out in the area surrounding an opencast quarry located in Vallerano near Rome. The plant extracts leucite basalt, that is a dark grey volcanic rock and has many domestic and industrial uses, the most common being its application as the raw material for railway track ballast. Figure 1 shows the site of the quarry and the surrounding residential area.
The mining site may be divided into two distinct areas: The face where the rock is mined with the help of explosives and the processing plant.
After blasting the fragments are extracted. The exploitation method is that of the single level splitting. Subsequently, the shot rock is taken by dump truck to the processing plant. Here it is crushed and sorted by various specific machines (crushers, grinding mills and screens) so that the required commercial granulometric sizes are obtained.
For the purposes of this paper, only the noise emissions generated by the processing plant are taken into consideration. The effects related to the use of explosives are mainly linked to ground vibrations which is beyond the scope of this paper.
From a methodological point of view the work was organized in successive phases. Firstly the physical model of noise propagation was employed starting with a survey of the data source. Subsequently a set of 50 measurements were taken inside the area of interest to be used as input for the geostatistical analysis.
The first step aimed at identifying and then classifying the fixed noise sources so as to determine the sound power level they each generate. For each of these machines a surface envelope was determined . Five measurements of the sound pressure were carried out with the machine under normal operating conditions. A Larson Davis model 824 sound level meter, which conforms to the requirements stated in D.M.16/03/1998, was used to measure the sound pressure level. Each measurement lasted 10 min.

Fig. 1. Quarry and surrounding residential area
On the basis of the measurements carried out, the sound power level for each machine is calculated as follows: Where: L w,f = Expresses the sound power level of the sound source in question at the single frequency [dB(A)] L p,f = The average sound pressure level, again for a single frequency, measured on the source surface envelope in real conditions [dB(A)] S and S 0 = Respectively the surface envelope of the source and the surface of reference set at 1 m 2 The average sound pressure level at the measurement point was determined while taking into account both the effects due to background noise and environmental effects (i.e., those related to the reflection of sound waves).
The average sound pressure level may be calculated thus: , With L′ average sound power level measured at the measurement points on the control surface, expressed in [dB(A)] and is defined as: With N number of measurements relative to each single source.
K 1 , correction factor, refers to the residual noise, i.e., the noise measured with the machine when off: This term includes all the other sources of noise other than those from the source being measured. Factor K 1 , expressed in dB, is given by: where, ∆L′ is the difference between the average sound pressure levels on the measured surface when the machine is on and when it is off. Lastly, K 2 , the environmental factor, is the term related to the increase in average sound level at the measured surface due to sound waves reflected by features of the environment. However, in open environments without reflective surfaces this term is null. The next step in determining the sound pressure level at the target point regards the application of the procedure proposed by the standard UNI EN ISO 9613 1-2, with appropriate calibrations for the present study.
The sound pressure level at the receptor, for each single frequency, is calculated by applying the following equation: , Where: L s,f = The predicted sound pressure level for the single frequency and at a certain distance from the sound emitting source [dB(A)] L w = The sound power level of the source (1) D I = The directivity index of the source K 0 is the shape factor relative to the emission characteristics of the source D s = The geometric divergence term defined as: where, d is the distance between source and target [m]; Lastly, the term inside the sum in (5), represents the set of factors which, for various reasons, contribute to the absorption of sound energy.
In particular: Where: D atm = The term related to absorption due to meteorological conditions [dB(A)] D g = The term related to absorption due to the ground [dB(A)] D b = That due to the presence of any barriers or obstacles [db(A)] The same assessment of environmental impact was carried out by means of geostatistical analysis. The basic concept of geostatistics is a probabilistic approach to modeling a physical mechanism. Each single value of a regionalized variable and above all its spatial variability structure, is assumed to be described efficiently in a stochastic framework. Data values are the results of many complex factors and can be reasonably viewed as possible outcomes of a random process. Therefore, from a probabilistic point of view, any spatial variable z(x), measured at a specific location, is the outcome of a random process that generated it from a random variable Z(x). At any point in the observed domain, the random variable may have different properties and thus generate different values. Hence these values are one possible outcome of the function that generates them following a certain regionalized probability density function. Considering that all the possible realizations are represented by the same random function, this function should be related to some spatial law. Such a law is represented by the cumulative distribution functions of n random variables.
Moreover, if the distribution function can be reasonably represented by the first two moments (e.g., Gaussian-like distribution) there is no need for a stationarity hypothesis except for the first two moments of the random function.
More specifically, these conditions may be expressed as follows: Equation 8 implies the first moment is spatially invariant, whereas (9) expresses covariance and does not depend on position but only on h distance between two selected points.
The variogram function defined as: Is a useful tool to describe spatial variability from sampled data. Once it is computed, if its shape is stationary, then the hypotheses (8) and (9) may be assumed. These characteristics, which are related to the random function connected with the phenomenon, may be appreciated, at least at a certain scale, when a sill is present in the variogram itself. Thus the variogram is expressed with regard to covariance as: Under these conditions, inference of the experimental variogram is performed and based on a least square fit of the experimental values of semivariance for each lag. Such a fit can be made by well-known continuous functions such as authorized models. Thus spatial variability is modelled and estimation may finally be carried out. The passage from a discrete information description to a continuous description is performed with the following linear estimation: In which coefficients are chosen to minimize the mean square prediction error (14) defined in (13), under the hypothesis of unbiasedness (15) where the term in square parentheses represents the difference between the sampled value and the estimated one at the same location: The previous equations are rewritten in a more compact way in the Kriging system: The only other notable feature is related to the choice of the regionalized variable. This is especially important since the estimation method is based on a linear interpolation of values (12), that cannot be performed with a logarithmic function, sound pressure level, defined a: Therefore the whole procedure is developed considering the regionalized variable as sound pressure (p a ) and only once the estimation is carried out, the conversion from sound pressure to its corresponding level is carried out.

Results and Discussion
In this study five sound sources were identified in the quarrying area (see red circles in Fig. 2) and the relative sound power level data are shown in Table 1. As expected, each source required ten measurements to be carried out on its own surface envelope. Of these, five were taken under normal operating conditions (the plant under maximum load and the source machine turned on), whereas the other five were taken under the same conditions but with the source machine turned off. Then the physical model according with Equation 5, was implemented in 50 points positioned within the area of interest at the locations illustrated in Fig. 2. At these points a measuring campaign was carried out monitoring each point for 30 min. The objective was to characterize the noise from the plant, such that in the postmeasurement phase, noises which originated from sources which are not directly correlated with the plant (voices of passersby, traffic sounds etc.) might be eliminated using suitable filters.
The sampled points were chosen as nodes of a regular rectangular grid covering the area of interest (Fig. 2). The distance between each node was set at 20 m with some exceptions due to local obstacles. This spatial step seemed to guarantee that the sample was not too large to characterize the acoustic climate. (Alfaro Degan et al., 2015).
The measured values were compared with those obtained from applying the model. The results are shown in Fig. 3 and 4. The parameters of the comparison indicate that the correlation coefficient between the two series of data is equal to 0.966.
The geostatistical analysis of the data was carried out on the same 50 samples in Fig. 2. The measurements (which lasted 30 min each) were carried out in order to fully characterize the temporal evolution of the signal by monitoring sound pressure level.
A systematic sampling strategy was used, defining a regular grid (the whole domain was 130×170 m) and sampling each node.
The first step was a structural analysis of the sampled data in order to identify the behaviour of the variable and any anisotropy.
Because of the importance of variographic computational parameters, polar diagrams were studied. On considering the characteristics of the sampling grid, no anisotropy was found in the available data. Furthermore, since the estimation should be calculated isotropically and the area should be related to the scale of the phenomenon, the computation domain was reduced to 20 lags.
The variogram, shown in Fig. 5, revealed a well defined spatial structure which was obtained by reducing the maximum lag distance to 120 m and thus following a quasi-stationary approach. The NW-SE direction was considered. Moreover, in order to represent the structural variability of the phenomenon, a fitting process of a mathematical function on the data variogram was performed.
This approach was developed with the continuous function (continuous line) which constituted a spherical model and a nugget effect.
Once the analytical variogram was defined, the estimation itself was carried out. Results are obtained by means of an Ordinary Kriging system (16) based on the nodes of a regular grid covering the entire domain.
A cross validation between the sampled values and the estimated values at each of the 50 available locations is shown in Fig. 6.
However, on analysing the system (16), it should be noted that at a point in which the regionalized variable is available and known, the Kriging system returns exactly the same value. In order to evaluate the quality of the results obtained, the 'leave one out' cross validation method is employed. This consists in considering the value of the variable at the point being tested as being unknown and estimating its value from the values of the remaining points in the set (without utilising the value of the point in question). The results are shown in Fig. 6 and 7.
Then, in order to analyse the quality of the results obtained, two parameters were considered: The correlation coefficient and the standard deviation of the differences between homologous values.
The data shown in Table 2 indicate a strong correlation in both cases and in both cases the standard deviation of the differences between the measured and calculated values was less than 2% i.e., 1.07 [dB(A)].    Then, in order to show the results geographically, two maps were plotted indicating the isobels from the data estimated from the Ordinary Kriging (Fig. 8) and the isobels from the output data of the ISO 9613 model (Fig. 9) respectively.

Conclusion
On analyzing the results it is apparent that with regard to a similar correlation coefficient for both methods, the ISO 9613 prediction had a greater standard deviation of differences and in at least two points displayed significant disparity due to local topological features such as abrupt changes in elevation or localized barriers (Fig. 3). The scatter plot in Fig. 5 does not show similar effects owing to the geostatistical method used. In general, with regard to the physical method, it may be said that implementing such a technique, requires a number of measurements which is dependent on the number of active sound sources. In this case 50 measurements were needed to characterize the 5 available sound sources and to define the background noise.
Clearly, as the number of active sound sources at the site increases, so the computational complexity and the burden of the measurement campaign become greater. On the other hand a geostatistical approach involves the constraint regarding the characterization of the spatial variability of the variable.
In the present study, with respect to the 50 measurements made on a regular grid of approximately 20 m each side, a structure of defined variability emerged which allowed quasi-stationary modelling of the phenomenon. However, it is not possible a priori to determine a minimum number of readings required for the implementation of a Kriging method.
It may be said that on the scale in question (over a distance of the order of hundreds of metres) and in reference to the data presented here, the geostatistical approach was successful: The number of measurements carried out being equal, acoustic climate values were obtained which fitted the readings taken in the field. It is therefore possible to conclude that in the case of industrial plants with inhomogeneous and variegated noise emissions (which is usually the case with extractive sites) geostatistical modelling may well be an alternative approach which is not only viable but also convenient.

Author's Contributions
Dario Lippiello: Designed the research plan, participated in all experiments, coordinated the data analysis and contributed to the writing of the paper.
Guido Alfaro Degan: Organized the study, participated in all experimets and contributed to the writing of the paper.
Mario Pinzari: Designed the research plan with Lippiello and critically reviewed the manuscript.

Ethics
All of the other authors have read and approved the manuscript and there are no ethical issues involved.