Urban air pollution have a direct impact on public health. Ultra ne particles (UFPs) are ubiquitous in urban environments, but their distribution are highly variable. In this paper, we take data from mobile deployments in Zurich collected over one year with over 25 million measurements to build a high-resolution map estimating the UFP distribution. More speci cally, we propose a new approach using a Gaussian Process (GP) to estimate the distribution of UFPs in the city of Zurich. We evaluate the prediction estimations against results derived from standard General Additive Models in Land Use Regression, and show that our method produces a good estimation for mapping the spatial distribution of UFPs in many timescales.
Tra c junctions, industrial installations and urban canyons all contribute to the high spatial and temporal variability of air pollution in urban areas. Small-scale spatial distribution of ambient air pollution have traditionally been studied with Land-Use-Regression (LUR) summarised in Hoek et al. (2008). It uses land-use and tra c characteristics of a particular grid region as explanatory variables to learn to estimate pollution concentrations under a Generalized Additive Model (GAM). The learnt model is then used to predict pollution levels for all locations with the available land-use information.
In this paper, we propose a novel approach of estimating urban ultra ne particle levels across di erent temporal aggregates from measurements collected from the trams. Similar to standard models in land use regression, it estimates the pollution levels within di erent grid-cells in the urban environment from a set of land-use features. Our model is based on constructing a Gaussian Process described in Rasmussen and Williams (2006), with additional consideration to spatial features in the covariance matrix. Following the practice in previous work of Hasenfratz et al. (2014), we evaluate the models (GAM, pure land use and mixed spatial-land-use) using standard random 10-fold cross validation.
The outline of this paper is as follows: we begin with a summary of the background to the paper: the data, the traditional models used in land-use regression, and introducing Gaussian Process Regression. We then introduce a new approach for estimating UFP levels, and evaluate it against the previous approach over a benchmark dataset. 2 2.1
The Aggregate Datasets
The data were selected from UFP measurements collected on Zurich trams between April 2012 and March 2013 as
part of the OpenSense project and the sensing methodology is described in Li et al. (2012b) an
Several preprocessing steps were used before the data were prepared for the model, including removal of
measurements within the indoor tram depot, measurements with bad GPS data, and measurements with
extraordinary high levels >100'000 particles per cm3. These steps were described in
We then aggregate the data within the di erent grid cells according to the di erent time windows, such
as yearly, seasonal, monthly, biweekly, weekly, daily and half-daily. This is done to understand the trade-o
between long and short term aggregate data. In order to evaluate and compare our results to previous work, we
followed the convention of selecting only the 200 grid cells with the highest measurements count for the purpose
of modelling and vali
Land Use Regression In literature, land-use regression models are used to assess intra-urban air pollution distributions, and a comprehensive review of these techniques can be found in Hoek et al. (2008). They typically combine monitoring of air pollution at 20-100 locations spread over the study area, and develop a model using predictor variables obtained through geographic information systems (GIS). The predictor variable generally include some tra c information, population density, designated land use and features of the landscape such as attitude and slope. Due to the cost of deployments, studies usually last 1-2 weeks in duration.
For particulate matter such as PM2:5 and PM10 and UFPs, Generalized Additive Models (GAMs) have been
used in land use regression to study their spatial and temporal variability. It typically use the following equation
to model the relationship between the pollution level p and a set of explanatory variables A1; : : : An.
ln(p) = a + s1(A1) + s2(A2) +
+ sn(An) +
(1)
where a is known as the intercept, the error term, and s1 : : : sn are typically smooth regression splines with
an upper limit of 3 on the degree of freedom. In this paper, we use the GAM
Gaussian Process Regression
Also known as Kriging, Gaussian process regression (GPR) has been extensively used for decades in Geostatistics
to model various spatial phenomena such as soil concentrations, weather-related events, etc., and in-depth
overviews can be found in Cressie and Cassie (1993) and Rasmussen and Williams (2006). Similar to other
non-parametric approaches, GPR does not require prior structural knowledge about the phenomenon. Indeed,
the idea is precisely that structure is directly inferred from the data. Furthermore, GPR outputs statistical
predictions and thus represents an adequate candidate to model phenomena that are inherently noisy and which
one can only observe through noisy instruments. Recently it has been successfully applied in many machine
learning tasks such as bioinformatics in Chu et al. (2005), sensor calibration in Monroy et al. (2012) and
crowdsourcing Venanzi et al. (2013) It still represents a very active ongoing research area as seen in e.g. Bonilla et al.
(2010), Cao et al. (2013) an
A Gaussian Process (GP) is used to model a phenomenon that takes place in a certain input space X Rd. We formally write f (x) where x 2 X the function that models the phenomenon. The general idea is to assume that the function f (x) is a speci c realization of a prior Gaussian Process GP, which is the generalization of a multivariate normal distribution to an in nity of random variables, that is to say a distribution over whole functions. A GP is fully de ned by its mean function m(x) and its covariance function k(x; x0) (also called kernel) that are the generalization of the mean vector, respectively the covariance matrix of a multivariate normal distribution.
Regression with a GP is typically performed as follows. In general, we can only make from the phenomenon noisy observations yi = f (xi) + i where the additive noise is also assumed to be Gaussian N (0; n2). By using the marginalization property of GPs and the additive nature of the noise we know the joint distribution of the observations y at locations X and the values f at test points X to be: y f
N
mm((XX)) ; k(Xk;(XX +;Xn)2I)
k(X; X )
k(X ; X )
Our rst GP model uses only land-use variables as features to generate predictions on the mean UFP
concentration measured by the sensors within the respective grid cells in the timeframe of the speci ed dataset. They
follow from the features use
m(x)) V[f (x )] = k(x ; x ) k(x ; X)(k(X; X) + n2I) 1k(X; x )
The main challenge is to create and choose prior mean and covariance functions that carry adequate assumptions about the phenomenon. We describe in detail how we derived such functions in the following section. k(xLU; x0LU) = f2 exp (xLU x0LU) where M = diag(`LU ) 2 (2) (3a) (3b) (4a) (4b) Swiss Federal Statistical O ce Canton of Zurich government
{ Average daily tra c volume OpenStreetMaps.org { Population density, industry density, building heights, heating type, terrain elevation, terrain slope { Main road type, distance to next major road, distance to major tra c signal
As we wanted to start with no particular a priori structural knowledge, only very simple mean functions were tried such as the trivial xed 0 function and a constant c. Deriving a suitable covariance function was, however, a bit more complex. Indeed, to be valid a covariance function must be positive de nite. It is common practice to start from well-known parametrized families of positive de nite functions and t the parameters (that in the scope of GPR are called hyperparameters) using the data. All the covariance functions that were tried are stationary that is to say every points of the space shows the exact same covariance structure with its own surroundings or more formally we have k(x; x0) = k(x x0). Stationary covariance functions such as squared exponential, and various avours of the Matern class were tried, each one carrying di erent assumption about the smoothness of the process. Finally from preliminary tests results, a constant mean function and a squared exponential covariance function were selected. We note that the chosen covariance function was the one that carried the strongest smoothness assumptions. The prior GP is thus de ned by:
m(xLU) = c The f2 is the magnitude hyperparameter, and the `LU are the length-scale hyperparameters that determine the relevance of some or other land-use variables. To learn the values of all the hyperparameters = (c; f2; `LU ; n2) one can either use optimization or sampling techniques. In our case, we used the standard approach that consists in optimizing the log marginal likelihood: Even though the explanation of the phenomenon given by the land-use variables may already be quite good, it is very likely that part of it still elude us because of some contributions to the phenomenon that are badly or not at all re ected in the variables. To address this matter, we tried to incorporate geographical informations into the model with the hope that such missed contribution will at least be partly explained locally.
The problem with parametric models such as GAM is that we cannot easily add geographic informations into the model in a sensible way. For example if we naively add the longitude and latitude as covariates, we would be making very strong assumptions rather unrealistic.
However, with GPR (and this is why it has been extensively used in Geostatistics) it is natural to include such informations in the reasoning. This is done by including a consideration for geographical distance in the covariance function. We call our second model a mixed spatial-land-use model, which is a variant of the rst one in which we added a term in the covariance structure. We also tried di erent isotropic kernels to be this additional term. From the preliminary experiments the following covariance function was selected: k( xxLSU ; xx0L0SU ) = f2LU exp It is worth noting that it is the exponential function, the less smooth of the considered covariance functions, that was chosen to be the additional term in function of the geographical distance. The values of the hyperparameters = (c; f2LU ; f2S ; `S ; `LU ; n2) were once again xed using marginal likelihood maximization. 4
We implemented our own Java framework to perform GPR. However, the conjugate gradient optimizer, used
to maximize the log marginal likelihood, was taken from the Matlab toolbox GPML v.2 (see Rasmussen and
Nickisch (2010)) and translated in Java. Most of linear algebra operations were carried out using EJML1 library.
The experiments were conducted on a server with 64 AMD Opteron processing cores and 96 GB of RAM. In the
experiments, we compared the following three di erent type of models on the UFP datasets described earlier.
1. GAM A General Additive Mo
From the benchmarking data supplie
Fig. 3 shows the satellite image of the urban area covered in the deployment, the output of the pollution map
for the season of summer in 2012, and the comparison of the prediction against ground truth of the same season
under random 10-fold cross validation. Fig.4 shows the scatter plots of model predictions against ground truth
1http://code.google.com/e cient-java-matric-library/
across all time scales, where all predictions of the same time scale are located on the same plot. It is worthy
to note that similar to the previous model presente
RMSE First we compare the Root Mean Square Error (RMSE) of the predictions derived from the models under random 10 fold cross validation (Fig. 5). It is a standard metric of predictive power for measuring the accuracy of prediction models. It is obtained by:
RM SE = s
PiN=1(pi
N gi)2 (6) where pi denotes the ith prediction, gi the ground truth of the ith prediction, and N the total number of predictions. In Fig. 5, the plot on the left displays the overall mean of the RMSE, while the box-plot on the right displays the minimum, lower quartile, median, upper quartile and maximum of the average RMSE of the whole 10-fold validation tests on all the datasets of the same time scale. The yearly data came from a single dataset, thus it is presented as a single value. It shows that as expected, the higher temporal resolution leads to higher uncertainty in the prediction, the GP models outperforms GAM across all temporal resolutions, and the mixed spatial-land-use model produced less error than the land-use only model. 3./]trcm 6000 a p [ ESM 4000 R 2000
0 2R 1 0.8 0.6 0.4 0.2 0 F2AC 0.9 1.2 1.1 1 0.8 0.7 0.6
GP_LU GP_LUXY GAM yearly seasonal monthly biweekly weekly daily halfdaily yearly seasonal monthly biweekly weekly daily halfdaily 2R 4.2
R2 Score
We then compare the R2 coe cient, also known as the coe cient of determination of the model predictions
(Fig. 6). It indicates how well the observed outcomes are replicated by the model predictions as the proportional
variation of outcomes explained by the model. Its formula is given by:
where pi denotes the ith prediction, gi the ground truth of the ith prediction, and g is the mean of the ground
truth. In Fig. 6 observe that the variance of the R2 score also increases as the the time scale shrinks across all
models. We see that the results of GP models in general have a higher R2 than GAMs, and introducing the
spatial covariance in the GP model also improves the R2 score across all time scales.
Finally, we compare the F AC2 score of the model predictions (Fig. 7). It measures the fraction of data points
that lie inside the factor of two area. It is a robust measure of prediction as it is not overly in uenced by high
and low outliers. It is derived by:
We implemented two schemes based on Gaussian Process for estimating mean UFP concentrations in urban areas
of Zurich, Switzerland. We show that they provide an alternative to GAM approaches in land-use regression,
and there is a general trade o between the length of the time scale and the quality of the model predictions. We
also show that across the timescales the proposed GP models presents an improvement on the current state of
the art. The resulting maps may be useful for application such as assessing population exposure to air pollutants
similar to that of Carroll et al. (1997), uncover areas of high air pollution for persons with allergies, or evaluate
the trustworthiness of measurements contributed by a community of sensors as described in Li et al. (2012a) an
Possible future work includes moving away from a grid-based model to make use of urban spatial features described in Li et al. (2012b), developing models that handles di erent aspects of sensor reliability and measurement bias; detecting and ltering spurious measurements, and combining meteorological information and real time data to produce the best real-time estimations for individual exposure analysis and route planning. Our approach based on Gaussian Process Regression is very general, and it is interesting to see if it can be generalised to particulate dispersion outside urban environments to applications such as bush- re detection; and whether it can be applied to estimating other air-borne or water-borne pollutant dispersions.
Acknowledgements We thank our collaborators at ETHZ David Hasenfratz and Olga Saukh for supplying the benchmarking data and model from their previous work. This work is supported by OpenSense project funded by NanoTera.ch, and the ARC Discovery Project (DP120103758) \Arti cial Intelligence Meets Sensor Networks". K. Aberer, M. Hauswirth, and A. Salehi. A middleware for fast and exible sensor network deployment. In
VLDB, 2006.
K. Aberer, S. Sathe, D. Chakraborty, A. Martinoli, G. Barrenetxea, B. Faltings, and L. Thiele. OpenSense:
Open community driven sensing of environment. In ACM IWGS, 2010.
Edwin V Bonilla, Shengbo Guo, and Scott Sanner. Gaussian process preference elicitation. In NIPS, pages 262{270, 2010.