Journal of Environmental Quality 31:121-128 (2002)
© 2002 American Society of Agronomy, Crop Science Society of America, and Soil Science Society of America
TECHNICAL REPORT
Atmospheric Pollutants and Trace Gases
Optimization of a Monitoring Network for Sulfur Dioxide
E. P. J. Boer*,a,
A. L. M. Dekkersb and
A. Steina
a Wageningen University & Research Centre, Mathematical and Statistical Methods Group, Dreyenlaan 4, 6703 HA,Wageningen, the Netherlands
b RIVM, National Institute of Public Health and the Environment, P.O. Box 1, 3720 BA, Bilthoven, the Netherlands
* Corresponding author (E.P.J.Boer{at}ato.dlo.nl)
Received for publication February 5, 2001.
 |
ABSTRACT
|
|---|
In this study, we develop and apply a methodology to reduce an existing monitoring network to find an optimal configuration of a smaller network. We use a criterion based on locally weighted regression with two different weight functions. The methodology is applied to the Dutch national SO2 network and offers the possibility to include different politically relevant options in the model by weight criteria. Because full enumeration of all monitoring networks is impossible, a combinatorial search algorithm is applied to find a (sub)optimal solution. If different (optimal) monitoring networks result from different criteria, then the best can be selected.
|
INTRODUCTION
|
|---|
IN ENVIRONMENTAL monitoring programs, optimization of the monitoring network is often an important issue. Most networks for air pollution were designed long ago, and it is sometimes necessary to re-examine a given network because the underlying physical processes and emissions may have changed so that collected data somehow do not answer the necessary questions. Several papers deal with optimization of monitoring networks (e.g., Caselton and Zidek, 1984; Warrick and Myers, 1987; Cressie et al., 1990; Pardo-Igúzquita, 1998; Van Groenigen and Stein, 1998; Fedorov et al., 1999). The goal of optimizing an environmental monitoring network is, in many cases, related to the accuracy of maps and/or reduction in costs.
Adaptation of an existing monitoring network is often done under constraints of authorities and sometimes based on expert judgment without any formal mathematical criteria. An example of a scientific criterion is the maximization of the minimum distance between monitoring stations, so that the stations are as evenly spread as possible over the area of interest (Müller, 1998). Other criteria are based on geostatistics such as minimization of the kriging variance (Cressie et al., 1990). If little is known about the structure of the stochastic process that underlies the monitoring data, locally weighted regression (Cleveland, 1979) is an appropriate interpolation technique. This flexible method allows the characterization of trends by using simple local models. The variance of estimates resulting from this method can be used to formulate a criterion for optimization.
The aim of this study is to further explore the possibilities of locally weighted regression for the optimization of a monitoring network (Müller, 1995; Fedorov et al., 1999). We show that it leads to a flexible criterion for adaptation of an existing monitoring network. This flexibility is expressed in the application of this interpolation technique and by incorporating different design criteria, which can easily include external (policy) constraints. The criterion is applied to a case study in the Netherlands, where the number of SO2 monitoring stations has to be decreased by about 60%. The combinatorial optimization problem of selecting the optimal monitoring network is solved by search algorithms.
|
LOCALLY WEIGHTED REGRESSION
|
|---|
Model Formulation
Let {x1,x2,...,xn} be locations, xi = 
, of n monitoring stations in a region D of the monitoring network 
n, with corresponding observations {y(x1),y(x2),...,y(xn)} on a certain point of time. The observations are modeled by:
 | [1] |
where 
(x) for x
D is a smooth function and 
i are independent, identically distributed, zero-mean observation errors. A flexible collection of functions consists of those functions that can be approximated locally by a polynomial expression. Consider a location x*D, and let L (x,ß(x*)) be the local approximation of (x), where ß(x*) is a vector with the parameters of a locally spatial trend around x* = 
. The smaller the distance between x and x* the better the approximation.
For local approximation, we consider the following polynomial expression of order p:
 | [2] |
which can be considered as a two-dimensional Taylor expansion, where ri is the remainder term of the approximation.
For an isotropic field, Eq. [2] can be simplified to:
 | [3] |
where h = ||x - x*||. For every estimation location x*, the vector ß(x*) has to be estimated. Weighted least squares is applied with decreasing weights as the distance between xi from x* increases. Let an estimator 
(x*) be defined as:
 | [4] |
where 
is a weight function depending on observation locations xi and estimation location x*. The use of a weight function corresponding to the model assumption that points close to x* plays a larger role in the determination of 
L (xi,ß(x*)) than points further away.
In principle, many weight functions can be considered. In this study, the choice is restricted to two weight functions mentioned by Müller (1998). The first is the so-called tricube weight function, defined as:
 | [5] |
where hf is a smoothing parameter, which determines the neighborhood of an estimation location x*. The McLain function is the second choice:
 | [6] |
where weights are only determined within a fixed neighborhood with range hf. Outside this neighborhood m 
= 0.
Given the polynomial expression, the weight function, and the smoothing parameter hf, the estimation of y(x*) is equal to:
 | [7] |
as can be seen from Eq. [2] and [3].
The smoothing parameter hf is estimated by calculating cross validation values for a range of smoothing parameter values. The optimal smoothing parameter is chosen as that value for hf that minimizes the following expression:
 | [8] |
where 
(xi)(-i) is the local fit at xi, without using y(xi) for the estimation, the "leaving-out-one method." The influence of the choice of the weight function will be shown in Results and Discussion, below.
Optimal Design for Locally Weighted Regression
Estimation variances of locally weighted regression parameters are used as a basis for a criterion for optimizing a monitoring network. These estimation variances are estimated on a set of locations where estimations are required: 
. Given a weight function and a value of the smoothing parameter hf, the classical theory of optimal design of experiments is applicable to every single estimation location. This allows formulation of a criterion for optimizing the monitoring network n = 
.
Let FTj be the design matrix in the case of an isotropic field and let p = 1 for the order of Eq. [3]. Then:
 | [9] |
where hij = ||xi - x*j||. If 
is diag(), = 
and the so-called information matrix is M-1j
= FTj
Fj, which depends on the design n, then:
 | [10] |
Since:
 | [11] |
only the first element of the parameter vector is required for estimation of y
. The estimated variance of 0 (x*) then equals:
 | [12] |
where 
11 is the upper left element of the matrix M-1j
and 
2 is the residual variance.
A criterion for assessing the performance of a network design n can be derived from Eq. [12] by calculating a weighted sum over a set of q estimation locations 
:
 | [13] |
The value of 
(n) is used as criterion for obtaining an optimal monitoring network n. The inclusion of the weights wj allows some locations to be considered more important than others, because of external or political reasons (see Specification of Different Design Criteria, below). An optimal design n is the design that minimizes the value of (n) in Eq. [13].
Case Study
In 1993, the RIVM (National Institute of Public Health and the Environment in the Netherlands) measured SO2 at 74 measuring stations of the Dutch National Air Quality Monitoring Network (Doesburg et al., 1994). Figure 1
shows the original monitoring stations with their annual average concentrations (µg m-3). To reduce the expense of data collection this network was reduced to 29 stations, a reduction of 60%. This reduction is feasible because SO2 concentrations are decreasing and the political pressure to maintain an expensive monitoring network is decreasing. In addition, deterministic models are available that allow the reliable calculation of SO2 concentrations (Bleeker and Den Hartog, 1995). However, a total abandoning of the entire network is not possible because of national and European regulations.
View larger version (31K):
[in this window]
[in a new window]
|
Fig. 1. Monitoring network of 74 stations in the Netherlands in 1993 coded with the annual average SO2 concentration (µg m-3).
|
|
Figure 1 shows that both the annual average SO2 concentrations and spatial variability of these values (i.e., relatively large differences in concentrations on short distance) are higher in the southwestern Netherlands compared with the northern Netherlands. Annual mean concentrations at 74 stations have an average of 10.2 µg m-3, whereas the minimum and maximum values are 4.2 µg m-3 and 28.1 µg m-3, respectively, and the variance is 26.6 µg2 m-6. A geostatistical approach toward optimization (Van Groenigen, 1999) was not an attractive option, because of difficulties in modeling the trend and because of the high spatial variability in the southwestern compared with the northern Netherlands. Therefore, the trend surface is estimated by a nonparametric regression technique: locally weighted regression, which was described in Model Formulation, above.
Specification of Different Design Criteria
The formulated criterion for the optimization of a monitoring network (Eq. [13]) is used for both weight functions in Eq. [5] and [6] to reduce the number of monitoring stations from 74 to 29. Given a particular choice for a weight function, the smoothing parameter hf has to be determined. By a reduction, however, the optimal smoothing parameter increases since the number of monitoring stations decreases. A small monitoring network leads to high variability of estimated values of the smoothing parameter. It is possible, however, to determine optimal smoothing parameters for several monitoring designs given the desired number of stations n. The average of the smoothing parameters thus obtained can be considered as an optimal smoothing parameter for a monitoring network of 29 stations.
Also, the number and coordinates of estimation locations have to be determined with corresponding weights (wj). Three sets of weights are used in this paper in order to accommodate external information or political regulations into Eq. [13]. The following criteria are considered:
- Criterion I puts equal weights at every estimation location.
- Criterion II puts weights according to the population density in the different provinces. It is an example of including political regulations into the optimization of a monitoring network.
- Criterion III puts weights in such a manner that it reflects the differences in spatial variability, that is, the greater variability in SO2 values in the southwestern compared with the northern Netherlands on short distances. The idea behind these weights is that in the north fewer monitoring stations are required than in the southwest because of the relatively constant values of the annual average SO2 concentrations. The weights are based on mean residual values of weighted local regression at the measurement locations in a neighborhood with range hf for every grid location. In this case, for p = 1, the isotropic, tricube weight function with hf = 93.
- Figure 2
shows estimation locations with corresponding weights for Criteria II and III.
View larger version (30K):
[in this window]
[in a new window]
|
Fig. 2. Estimation locations with corresponding weights based on population density at each location in number of inhabitants per km2 (Criterion II, left). Estimation locations with weights based on residuals of locally weighted regression (Criterion III, right).
|
|
Optimal Reduction of a Monitoring Network
The reduction of an existing monitoring network is a combinatorial optimization problem (Ko et al., 1995). For the case study we have to select 29 stations from a full set of 74 possible monitoring stations. This combinatorial optimization problem yields a large number of possible combinations (3 x 1020). It is impossible to enumerate these in full to find the best solution in terms of Eq. [13] and the different design criteria. Therefore, we developed a sequential search algorithm to find a (sub)optimal solution (see Sequential Search Algorithm, below). In Search Algorithm with Combinatorial Analysis, below, this algorithm is improved by using combinatorial solutions of subproblems.
Sequential Search Algorithm
A combination of a drop and an add algorithm is applied as a search algorithm. The drop algorithm sequentially removes stations from a set of monitoring stations k until n monitoring stations are left. Mostly the algorithm will be started at k = N, where N is the number of stations in the original monitoring network. The set of removed points is defined as +N-k, which is used later on in the drop–add algorithm. The drop algorithm is described as follows (Rasch et al., 1997):
Drop algorithm
- Let k =
and +N-k = 
.
- Determine * = min
over t = 1,...,k.
- Drop the point xt corresponding to * from the design k; add the point to +N-k; k := k - 1; renumber set k.
- If k = n, STOP; otherwise go to 2.
The add algorithm can be formulated in a manner analogous to the drop algorithm (Rasch et al., 1997). The idea of the drop–add algorithm is that some points can be exchanged after running the drop algorithm. The number of additional points added to, or dropped from, a design is called m. The drop–add algorithm consists of three steps. First, an (n - m)-point design is obtained with the drop algorithm, followed by an addition of 2 x m points and finally deletion of m points, so that an n-point monitoring design is obtained. This is repeated until no improvements are found.
Search Algorithm with Combinatorial Analysis
Although a full enumeration of all combinations is impossible, smaller combinatorial subproblems can be solved. Sequential search algorithms add or drop one point at each iteration.
It may be worthwhile to consider dropping combinations of points. Rasch et al. (1997) describe a branch-and-bound algorithm that solves combinatorial problems for optimal designs in regression analysis. With some adaptations, the branch-and-bound algorithm can be applied to the optimization problem discussed in this paper. The branch-and-bound algorithm is part of the last step of the drop–add algorithm, which is the final dropping of m points.
|
RESULTS AND DISCUSSION
|
|---|
Maps of Interpolated Sulfur Dioxide Values with Locally Weighted Regression
As pointed out in Model Formulation, above, three choices have to be made in order to apply locally weighted regression. For the order p of the polynomial expression, p = 1, both isotropic and anisotropic are decided upon first, being a compromise between computational ease and flexibility (Müller, 1998). A higher order of the polynomial expression did not improve the prediction accuracy. The McLain and the tricube weight function—Eq. [5] and [6]—were used as a weight function. Given a weight function, the unknown smoothing parameter hf is estimated by cross validation. Figure 3
shows the cross validation results for p = 1, both isotropic and anisotropic for both weight functions and these weight functions under the optimal smoothing parameter. A simplification of the problem to an isotropic approach results in higher values of the cross validation. In the anisotropic case, a small range of a neighborhood hf is preferred for both weight functions. However, hf cannot be chosen too small because at least three measurement points have to be available in the neighborhood. Therefore, a range of a neighborhood of 70 km (hf = 70) is chosen for both weight functions. Furthermore, Fig. 3 shows that cross validation values for the McLain weight function are smaller for every hf (separately for isotropic and anisotropic). This indicates that locally weighted regression with the McLain weight function produces estimations with a higher prediction accuracy. However, caution is required because the results are only based on observations made at 74 monitoring stations. Maps of interpolated SO2 values are shown in Fig. 4 .
View larger version (16K):
[in this window]
[in a new window]
|
Fig. 3. Plots of cross validation values (N = 74) for a range of smoothing parameters for the tricube and the McLain weight function, both isotropic as anisotropic (top). The weight functions with the optimized smoothing parameters (bottom).
|
|
View larger version (42K):
[in this window]
[in a new window]
|
Fig. 4. Maps of interpolated SO2 values (N = 74) of the annual average concentration of SO2 by locally weighted regression, hf = 70. Tricube weight function (top) and McLain weight function (bottom).
|
|
The maps of interpolated SO2 values in Fig. 4 show a clear difference. The surface of interpolated SO2 values using the tricube weight function is smoother than the map obtained by using the McLain weight function, as a result of different shapes of the weight functions (Fig. 3, bottom). The McLain weight is almost an exact interpolator, because the point closest to x* gets a large weight compared with the weights of points further away.
Comparison of Search Algorithms
In Optimal Reduction of a Monitoring Network, above, three sequentially search algorithms are introduced: the drop algorithm, the drop–add algorithm, and the drop–add algorithm with combinatorial analysis. In this paragraph algorithms will be compared; because of computational reasons the calculations are restricted to the isotropic case. A smoothing parameter of 175 (km) for the tricube weight function and 183 (km) for the McLain weight function is chosen as a test case.
Figure 5
shows how the values of (n) increase when monitoring stations are dropped sequentially from an existing monitoring network. The values of (n) are calculated for the tricube weight function with a smoothing parameter of 175 km with equal weights at the q estimation locations. The smaller the monitoring network becomes, the larger the influence on the value of (n) of dropping one station from the monitoring design.
The drop algorithm is further refined to the drop–add algorithm and the drop–add algorithm with combinatorial analysis. The question that we would like to discuss in this paragraph is how far designs found by the sequential drop–add algorithm (A) differ from those found by the drop–add algorithm with combinatorial analysis (B). Therefore, we introduce efficiency of a monitoring design by dividing the value of (n) of the results of Algorithm B by (n) of the designs found by maximizing the minimum distance (maxmin distance design). This is a very simple way of calculating an optimal design, which is comparable with common practice. Table 1 presents these efficiencies. Note that only values of (n) of designs calculated with the same weight function and criterion (I, II, or III) can be mutually compared.
Table 1 shows that there is only a small difference between the drop–add algorithm with combinatorial analysis and the sequential drop–add algorithm. For the McLain weight function with Criterion III, the value of (n) was even slightly higher with Algorithm B (0.01%). The maxmin distance designs have considerably higher values of (n). The computation time for Algorithm B is approximately 80 min, about 200 times as much as for Algorithm A.
Optimized Monitoring Networks
As already pointed out in Specification of Different Design Criteria, a smoothing parameter for a given weight function has to be chosen before starting the optimization. The optimal values for the smoothing parameters found in Maps of Interpolated Sulfur Dioxide Values with Locally Weighted Regression, above, were obtained using the full network of 74 stations. The values of the smoothing parameters will be larger for reduced monitoring networks. To obtain an optimal smoothing parameter for a monitoring network of 29 stations, the 74 stations are devided into three clusters and a random selection of 29 stations (in total) is chosen from these clusters. This is followed by determination of the optimal smoothing parameter (i.e., hf with the lowest value of the cross validation). This is done 500 times, and the average of the optimal smoothing parameters is considered as the optimal hf for a monitoring network of 29 stations for all three criteria. A smoothing parameter of 120 (km) for the tricube weight function and 136 (km) for the McLain weight function is obtained by this procedure.
Results of the drop–add algorithm (Sequential Search Algorithm, above) when n = 29 and m = 6 for the three different design criteria are presented in Fig. 6
. A comparison of results for two weight functions shows that the McLain weight function tends to spread the monitoring stations more over the Netherlands than the tricube weight functions. Further, a clear influence of the different criteria (weights at estimation locations) can be seen in the configuration of monitoring stations, especially for the McLain weight function. Monitoring stations are moved from parts of the Netherlands with low weights to parts with high weights at estimation locations.
View larger version (32K):
[in this window]
[in a new window]
|
Fig. 6. Monitoring networks (black squares) obtained from the network of 74 monitoring stations by the drop–add algorithm  . Tricube weight function, p = 1, anisotropic, and hf = 120 (km), for the three design criteria (A–C) (see the Specification of Different Design Criteria section). Same for McLain weight function, hf = 136 (km) (D–F).
|
|
Choice of a weight function for local weighted regression has a considerable influence on the final optimal monitoring networks. The two weight functions considered in this paper have a different shape. The difference between weights put at points close to x* and points further away is small for the tricube weight function compared with the McLain weight function. The McLain weight function puts relatively high weights at points close to x* and weights diminish rapidly as the distance to x* gets larger. Therefore, the McLain weight function will tend to spread the monitoring stations as evenly as possible over the region.
The choice of the weight function is, in principle, arbitrary. An investigation with cross validation can help to choose the best concerning the estimation accuracy. The smoothness of the surface of SO2 values is also controlled by the weight function. If there is interest in the mean SO2 concentration over the Netherlands a more smooth result will be desirable (tricube weight function). If the maximum SO2 concentration is of interest, a weight function as the McLain weight function will be preferred.
In this study we used locally weighted regression for interpolation and optimizations of an existing monitoring network. Other interpolation techniques such as stratified kriging could be used, wherein the southwestern Netherlands can be considered as a separable stratum. However, optimizing a monitoring network for the entire country can be a problem, because there are not enough monitoring stations in a stratum to fit a semivariance function to use for kriging. Furthermore, it would be difficult to optimize a monitoring network near the boundary between two strata.
|
CONCLUSIONS
|
|---|
The reduction of a monitoring network is investigated in this study. Design of experiments for locally weighted regression is used to reduce the number of stations in an optimal way. Different criteria are formulated on the basis of different objectives for monitoring SO2, these are reflected in the weights assigned at estimation locations. We make the assumption throughout this study that is constant over the whole country. Possibly, if more data are available, then a locally estimated - (x) could be included into Eq. [13].
Two search algorithms are used: a sequential drop–add algorithm (A) and a drop–add algorithm with combinatorial analysis (B). Algorithm B requires more computation time than Algorithm A, with a slight improvement of efficiency (Table 1).
Different design criteria, distinguished by the choice of the weights at estimation locations, result in different optimal monitoring networks. A precise formulation of the monitoring objective is necessary to make sure that the optimized monitoring network is indeed adequate. Table 1 shows that the values of (n) of maxmin distance designs can be considerably greater than those of optimized designs. Although the objective of a monitoring network is difficult to formulate in practice, it should receive more attention.
|
ACKNOWLEDGMENTS
|
|---|
The authors would like to thank the Air Research Laboratory (LLO) of the RIVM for providing the data. They also thank A.C. van Eijnsbergen and E.M.T. Hendrix for giving helpful suggestions for this paper.
|
REFERENCES
|
|---|
- Bleeker, A., and P.R. Den Hartog. 1995. Rapportage besluiten luchtkwaliteit over het jaar 1993. Onderzoeksrapport VROM, Publikatiereeks lucht & energie. (In Dutch.) Staatuitgeverij, Den Haag, the Netherlands.
- Caselton, W.F., and J.V. Zidek. 1984. Optimal monitoring network designs. Stat. Prob. Lett. 2:223–227.
- Cleveland, W.S. 1979. Robust locally weighted regression and smoothing scatterplots. J. Am. Stat. Assoc. 74(368):829–836.[ISI]
- Cressie, N., C.A. Gotway, and M.O. Grondona. 1990. Spatial prediction from networks. Chemometr. Intellig. Lab. Syst. 7:251–271.
- Doesburg, M.J., E.C.M. Rentinck, and P. Swaan. 1994. Landelijk Meetnet Luchtkwaliteit. Meetresultaten 1993. Deel 1-4, Onderzoeksrapport RIVM, Rapport 722101010. (In Dutch.) RIVM, Bilthoven, the Netherlands.
- Fedorov, V.V. 1989. Kriging and other estimators of spatial field characteristics (with special reference to environmental studies). Atmos. Environ. 23(1):175–184.
- Fedorov, V.V., G. Montepiedra, and C.J. Nachtsheim. 1999. Design of experiments for locally weighted regression. J. Stat. Plan. Inference 81:363–382.
- Ko, C.W., J. Lee, and M. Queyranne. 1995. An exact algorithm for maximum entropie sampling. Oper. Res. 43(4):684–691.
- Müller, W.G. 1995. Optimal design for local fitting. J. Stat. Plan. Inference 55:389–397.
- Müller, W.G. 1998. Collection spatial data—Optimum design of experiments for random fiels. Physica-Verlag, Heidelberg.
- Pardo-Igúzquita, E. 1998. Optimal selection of number and location of rainfall gauges for areal rainfall estimation using geostatistics and simulated annealing. J. Hydrol. 210:206–220.
- Rasch, D.A.M.K., E.M.T. Hendrix, and E.P.J. Boer. 1997. Replication-free optimal designs in regression analysis. Comput. Stat. 12: 19–52.
- Van Groenigen, J.W. 1999. Constrained optimisation of spatial sampling, a geostatistical approach. Ph.D. diss. ITC Publ. Ser. 65. ITC, Wageningen, the Netherlands.
- Van Groenigen, J.W., and A. Stein. 1998. Constrained optimisation of spatial sampling using continuous simulated annealing. J. Environ. Qual. 27:1078–1086.[Abstract/Free Full Text]
- Warrick, A.W., and D.E. Myers. 1987. Optimization of sampling locations for variogram calculations. Water Resour. Res. 23(3):496–500.
TOP
ABSTRACT
INTRODUCTION
