Imprecise Uncertainty Modelling of Air Pollutant PM10

1. Introduction

Particulate matter (PM) refers to solid particles and liquid droplets found in air. Many manmade and natural sources produce PM directly, or produce pollutants that react in the atmosphere to form PM. The resultant solid and liquid particles come in a wide range of sizes, and particles that are 10 micrometers or less in diameter (PM₁₀) can be inhaled into and accumulate in the respiratory system and are believed to pose health risks (Environmental Protection Agency, 2010). Particulate matter is one of the six primary air pollutants the Environmental Protection Agency (EPA) regulates, due to exposure to high outdoor PM₁₀ concentrations causes increased disease and death (Environmental Protection Agency, 2010). Therefore, PM₁₀ concentrations, amongst many other air pollutants, are sampled and measured in various places in California, United States.

The general trend of PM air pollutant concentrations in the air in California are on the decrease, but it continues to be monitored and observed. The California standards for annual PM₁₀ concentrations is that the annual arithmetic mean is 20 µg/m³, and the national standard is 50 µg/m³ before 2006 (California Environmental Protection Agency Air Resources Board, 2010, Environmental Protection Agency, 2010). The State of California sets very high standards for their air quality, and air pollutants are carefully monitored.

However, in reality, it is too costly in terms of time, finance, and manpower to keep all the 213 sites to be monitoring and recording. In Fig. 1, a complete map of all 213 sample locations for PM₁₀ are shown. However, one must note that these sample sites are never all used at any given year, PM₁₀ samples are taken at different locations each year. At best, a maximum of 102 PM₁₀ samples are collected during some years, and at worst, 61 PM₁₀ samples are collected at that year. Therefore comparisons of PM₁₀ between years are difficult, due to missing data at sample sites. It is difficult to construct kriging maps in terms of actual observations annually since the air pollutants were measured in different locations each year although the site design originally planned was quite delicate statistically.

Each year, approximately 40% of the 213 sites were actually observed. We call a site that does not have a recorded PM₁₀ value as "missing value", and since there are no patterns so that serious problems would twist the kriging map constructions. In Fig. 2, this is clearly demonstrated. In 1989, there are 61 PM₁₀ samples collected (29% of 213 locations), and in 2000, there are 94 PM₁₀ samples collected (44%).

Figure 1.

Complete 213 Observational Sites in the California State

Figure 2.

PM₁₀ Samples Collected in California in 1989 (61 sites) and 2000 (94 sites)

The data scarcity brings in a series of (five) fundamental issues into the spatial-temporal modelling and prediction practices for California PM₁₀ data, namely:

1. The necessity to recognize the impreciseness in analyzing the spatial-temporal pattern in terms of California PM₁₀ records, which inevitably acts the solidness of a geo-statistical analysis;

2. Which theoretical foundations are appropriate for modelling impreciseness uncertainty;

3. How to fill up the "missing value" sites so that the "complete" records are available, which is either an original annual average from the original observations (40%) recorded on the site or a or predicted value by "neighbourhood sites" (60%), i.e., to facilitate spatial-temporal imprecise PM₁₀ value by interpolations and extrapolations;

4. How to estimate the parameters of uncertain processes (temporal patterns), particularly the rate of change parameter${\alpha }_i,i=1,2,\cdots ,213$;

5. Create annual kiging maps (19 maps) under spatially isotropy and stationarity assumptions so that the changes between annual maps can be analyzed by kriging map difference between 2007 and 1989 and kriging map of location rate of change.

These issues will be addressed in the remaining sections sequentially.

2. The necessity of modelling impreciseness in California PM₁₀ spatial-temporal analysis

Impreciseness is a fundamental and intrinsic feature in the PM₁₀ spatial-temporal modelling, due to the observational data shortage and incompleteness. Spatially, there are 213 sites involved, and temporally, PM₁₀ observations were collected from 1989 to 2007, over a 19-year period. During the 19-year period, there are only two sites (Site 2125 and Site 2804) having complete 19 year records. There are 16 sites having only 1 record (8%) and 70 sites having 10 or above records. To have a statistically significant time-series analysis, 50 data points are minimal requirement for each site, so classical time-series analysis (probabilistic analysis) cannot be performed. In order to have a quick overall evaluation of PM₁₀ records on each site, we borrow the statistical quality control idea here (Electric, 1956, Montgomery, 2001). But we do not carry on traditional 6-sigma rule, rather, classify the PM₁₀ records into four groups: 1-$\left(5,20\right]$, 2-$\left(20,35\right]$, 3-$\left(35,50\right]$, 4-$\left(50,160\right]$. These four-group limits in Table 1 reflect the national standard, (50 µg/m³) and California state standard (20 µg/m³) respectively. For example, 1 -$\left(5,20\right]$is for a location whose PM₁₀ fall in 5 to 20 (µg/m³).

County name	1-$\left(5,20\right]$	2-$\left(20,35\right]$	3-$\left(35,50\right]$	4-$\left(50,160\right]$	No. of Sites
Los Angeles			3	7	10
Kern	1	2	1	5	9
Riverside	3		3	3	9
San Diego			5	3	8
Imperial	3			4	7
Lake		3			3
Inyo	2	2	4		8

Table 1.

PM₁₀ Hazard level evaluation over selected 7 counties

One must be aware that the classification is not in absolute sense, rather, additional rules are adding (similar to quality control chart pattern analysis (Electric, 1956):

(1) if a single point, then, classify the site hazard level according to which group it falls in; (2) if a sequence of records, some of them, particularly early points may fall in higher (or lower) hazard level, but if last three points fall in a lower (or higher) hazard level, the later level would be chosen for the site.

The additional rule 1 can attribute to expert's knowledge confirmation, while the additional rule 2 can be regarded as an expert's decision based on trend pattern.

Figure 3.

The 7 sites from the selected 7 counties with original PM₁₀ data plots and the hazard level classifications

Fig. 3 shows the classifications of a seven sites from the selected 7 counties in Table 1, each county one site is picked up for illustration purpose. The red coloured plot means the hazard level 1$\left(5,20\right]$; the green coloured plot means the hazard level 2$\left(20,35\right]$; the purple coloured plot means the hazard level 3$\left(35,50\right]$; and the black coloured plot means the hazard level 4$\left(50,160\right]$.

It is evident that facing the impreciseness caused by incomplete data recording, one has to rely on expert's knowledge to compensate the inadequacy and accuracy in collected observational data. Impreciseness is referred to a term with a connotation specified by an uncertain measure or an uncertainty distribution for each of the actual or hypothetical members of an uncertainty population (i.e., collection of expert's knowledge). An uncertain process is a repeating process whose outcomes follow no describable deterministic pattern, but follow an uncertainty distribution, such that the uncertain measure of the occurrence of each outcome can be only approximated or calculated.

The uncertainty modelling without a measure specification will not have an rigorous mathematical foundations and consequently the modelling exercise is baseless and blindness. In other words, measure specification is the prerequisite to spatial-temporal data collection and analysis. For example, without Kolmogrov's (1950) three axioms of probability measure, randomness is not defined and thus statistical data analysis and inference has no foundation at all.

Definition 2.1: Impreciseness is an intrinsic property of a variable or an expert's knowledge being specified by an uncertain measure.

It is therefore inevitably to seek appropriate form of uncertainty theory to meet the impreciseness challenges. In the theoretical basket, interval uncertainty theory (Moore, 1966), fuzzy theory (Zadeh, 1965, 1978), grey theory (Deng, 1984), rough set theory (1982), upper and lower provisions (or expectations) (Walley, 1991), or Liu’s uncertainty theory (2007, 2010) may be chosen.

While imprecise probability theory (Utikin and Gurov, 1998) may be a typical answer to address the observational data inaccuracy and inadequacy. However the imprecise probability based spatial modelling requires too heavy assumptions. Just as Utikin and Gurov (2000) commented, “the probabilistic uncertainty model makes sense if the following three premises are satisfied: (i) an event is defined precisely; (ii) a large amount of statistical samples is available; (iii) probabilistic repetitiveness is embedded in the collected samples. This implies that the probabilistic assumption may be unreasonable in a wide scope of cases.” Guo et al. (2007) and Guo (2010) did attempt to address the spatial uncertainty from the fuzzy logic and later Liu's (2007) credibility theory view of point.

Nevertheless, Liu’s (2007, 2010) uncertainty theory is the only one built on an axiomatic uncertain measure foundation and fully justified with mathematical rigor. Therefore it is logical to engage Liu’s (2007, 2010, 2011) uncertainty theory for guiding us to understand the intrinsic character of imprecise uncertainty and facilitate an accurate mathematical definition of impreciseness in order to establish the foundations for uncertainty spatial modelling under imprecise uncertainty environments.

3. Uncertain measure and uncertain calculus foundations

Uncertainty theory was founded by Liu in 2007 and refined in 2010, 2011. Nowadays uncertainty theory has become a branch of mathematics.

A key concept in uncertainty theory is the uncertain measure, which is a set function defined on a sigma-algebra generated from a non-empty set. Formally, let $\Xi$ be a nonempty set (space), and $$\mathfrak{A}\left(\Xi \right)$$ the $\sigma$-algebra on$\Xi$. Each element, let us say,$A\subset \Xi$ ,$A\in \mathfrak{A}\left(\Xi \right)$ is called an uncertain event. A number denoted as$$ƛ\left\{A\right\}$$, $$0\le ƛ\left\{A\right\}\le 1$$, is assigned to event$A\in \mathfrak{A}\left(\Xi \right)$, which indicates the uncertain measuring grade with which event $A\in \mathfrak{A}\left(\Xi \right)$ occurs. The normal set function$$ƛ\left\{A\right\}$$satisfies following axioms given by Liu (2011):

Axiom 1: (Normality)$$ƛ\left\{\Xi \right\}=1$$.

Axiom 2: (Self-Duality) $ƛ\{\cdot \}$is self-dual, i.e., for any$A\in \mathfrak{A}\left(\Xi \right)$,$ƛ\left\{A\right\}+ƛ\left\{A^c\right\}=1$.

Axiom 3: ($\sigma \text{-}$Subadditivity) $$ƛ\left\{{\displaystyle \underset{i=1}{\overset{\infty }{\cup }}{A}_i}\right\}\le {\displaystyle \sum _{i=1}^{\infty }ƛ\left\{A_i\right\}}$$for any countable event sequence$\left\{A_i\right\}$.

Axiom 4: (Product Measure) Let $\left({\Xi }_k,{\mathfrak{A}}_{{\Xi }_k},{ƛ}_k\right)$be the $k^{th}$uncertain space,$k=1,2,\cdots ,n$. Then product uncertain measure $ƛ$on the product measurable space$\left(\Xi ,{\mathfrak{A}}_{\Xi }\right)$is defined by

$$ƛ={ƛ}_1\wedge {ƛ}_2\wedge \cdots \wedge {ƛ}_n=\underset{1\le k\le n}{\min }\left\{{ƛ}_k\right\}$$E1

where

$$\Xi ={\Xi }_1\times {\Xi }_2\times \cdots \times {\Xi }_n={\displaystyle \prod _{k=1}^n{\Xi }_k}$$E2

and

$${\mathfrak{A}}_{\Xi }={\mathfrak{A}}_{{\Xi }_1}\times {\mathfrak{A}}_{{\Xi }_2}\times \cdots \times {\mathfrak{A}}_{{\Xi }_n}={\displaystyle \prod _{k=1}^n{\mathfrak{A}}_{{\Xi }_k}}$$E3

That is, for each product uncertain event $\Lambda \in {\mathfrak{A}}_{\Xi }$ (i.e,$\Lambda ={\Lambda }_1\times {\Lambda }_2\times \cdots \times {\Lambda }_n\in {\mathfrak{A}}_{{\Xi }_1}\times {\mathfrak{A}}_{{\Xi }_2}\times \cdots \times {\mathfrak{A}}_{{\Xi }_n}={\mathfrak{A}}_{\Xi }$), the uncertain measure of the event $\Lambda$is

$$ƛ\left\{\Lambda \right\}=\{\begin{array}{c}\begin{array}{cc}\begin{array}{cc}\underset{A_1\times \cdots \times A_n\subset \Lambda }{\sup }& \underset{1\le k\le n}{\min }ƛ\left\{{\Lambda }_k\right\}\end{array}& \text{ if }\begin{array}{cc}\underset{A_1\times \cdots \times A_n\subset \Lambda }{\sup }& \underset{1\le k\le n}{\min }ƛ\left\{{\Lambda }_k\right\}\end{array}>0.5\end{array}\\ \begin{array}{cc}\begin{array}{cc}1-\underset{A_1\times \cdots \times A_n\subset {\Lambda }^c}{\sup }& \underset{1\le k\le n}{\min }ƛ\left\{{\Lambda }_k\right\}\end{array}& \text{if }\begin{array}{cc}\underset{A_1\times \cdots \times A_n\subset {\Lambda }^c}{\sup }& \underset{1\le k\le n}{\min }ƛ\left\{{\Lambda }_k\right\}>0.5\end{array}\end{array}\\ \begin{array}{cc}0.5& \text{ otherwise}\end{array}\end{array}$$E4

Definition 3.1: (Liu, 2007, 2010, 2011) A set function $$ƛ:\mathfrak{A}\left(\Xi \right)\to \left[0,1\right]$$ satisfies Axioms 1-3 is called an uncertain measure. The triple $\left(\Xi ,\mathfrak{A}\left(\Xi \right),ƛ\right)$ is called an uncertainty space.

Definition 3.2: (Liu, 2007, 2010, 2011) An uncertainty variable is a measurable function $\xi$from an uncertainty space $\left(\Xi ,\mathfrak{A}\left(\Xi \right),ƛ\right)$to the set of real numbers, i.e., for any Borel set B of real numbers, the set$\left\{\tau \in \Xi :\xi \left(\tau \right)\subset B\in \mathfrak{B}\left(ℝ\right)\right\}\in \mathfrak{A}\left(\Xi \right)$, i.e., the pre-image of B is an event.

Remark 3.3: Parallel to revelation of the connotation of randomness in geostatistics, impreciseness occupies an fundamental position in geospatial-temporal uncertainty statistical analysis. In California PM₁₀ spatial-temporal study, nearly 60% sites do not have "complete" temporal sequences so that in order to fill the "missing" observations, we have to engage expert's knowledge to pursue "complete sequences" (i.e., to have 19 PM₁₀ values at each individual site), which is inevitably imprecise and incomplete. Impreciseness is referred to a term here with an intrinsic property governed by an uncertainty measure or an uncertainty distribution for each of the actual or hypothetical members of an uncertainty population (i.e., collection of expert's knowledge). An uncertainty process is a repeating process whose outcomes follow no describable deterministic pattern, but follow an uncertainty distribution, such that the uncertain measure of the occurrence of each outcome can be only approximated or calculated.

Remark 3.4: Impreciseness exists in engineering, business and research practices due to measurement imperfections, or due to more fundamental reasons, such as insufficient available information,..., or due to a linguistic nature, because it is an unarguable fact that impreciseness exists intrinsically in expert’s knowledge on the real world.

Definition 3.5: Let$\xi$be a uncertainty quantity of impreciseness on an uncertainty measure space$\left(\Xi ,\mathfrak{A}\left(\Xi \right),ƛ\right)$. The uncertainty distribution of $\xi$ is

$${\Psi }_{\xi }\left(x\right)=ƛ\left\{\tau \in \Xi |\xi \left(\tau \right)\le x\right\}$$E5

An imprecise variable $\xi$ is an uncertainty variable and thus is a measurable mapping, i.e.,$\xi :\mathbb{D}\to ℝ,\mathbb{D}\subseteq ℝ$. An observation of an imprecise variable is a real number, (or more broadly, a symbol, or an interval, or a real-valued vector, a statement, etc), which is a representative of the population or equivalently of an uncertainty distribution${\Psi }_{\xi }(\cdot )$ under a given scheme comprising set and $\sigma$-algebra. The single value of a variable with impreciseness should not be understood as an isolated real number rather a representative or a realization from the uncertain population.

Definition 3.6: (Lipschitz condition) Let $f\left(x\right)$be a real-valued function,$f:ℝ\to ℝ$. If for any$x,y\in {ℝ}^n$, there exists a positive constant$M>0$, such that

$$\left|f\left(x\right)-f\left(y\right)\right|<M\left|x-y\right|$$E6

Definition 3.7: (Lipschitz continuity) Let $f:{ℝ}^m\to {ℝ}^m$

1. $\forall B\subset {ℝ}^m$$B$$f$$B$$\exists M>0$

$$\Vert f\left(x\right)-f\left(y\right)\Vert <M\Vert x-y\Vert ,\forall x,y\in B$$E7

where $\Vert \cdot \Vert$is some metric (for example, Euclidean distance in${ℝ}^m$), such

$$d\left(f\left(x\right),f\left(y\right)\right)<Md\left(x,y\right),\forall x,y\in B$$E8

1. for each$z\in {ℝ}^m$, $f$is Lipschitz continuous locally on the open ball$B$of center $z$radius $M>0$ such

$$B_M\left(z\right)=\left\{y\in {ℝ}^m|d\left(y,z\right)<M\right\}$$E9

1. $f$${ℝ}^m$

Remark 3.8: For continuity requirements, Lipschitz continuous function is stronger than that of the continuous function in Newton calculus but it is weaker than the differentiable function in Newton differentiability sense. In other words, Lipschitz-continuity does not warrant the first -order differentiability everywhere but it does mean nowhere differentiability. Lipschitz-continuity does not guarantee the existence of the first-order derivative everywhere, however, if exists somewhere, the value of the derivative is bounded since

$$\frac{\left|f\left(x\right)-f\left(y\right)\right|}{\left|x-y\right|}<M$$E10

by recalling the definition of the Newton derivative

$$\underset{y\to x}{\lim }\frac{f\left(x\right)-f\left(y\right)}{x-y}=f\text{'}\left(x\right)$$E11

Similar to the concept of stochastic process in probability theory, an uncertain process $\left\{{\xi }_t,t\ge 0\right\}$is a family of uncertainty variables indexed by $t$and taking values in the state space$\mathbb{S}\subseteq ℝ$.

Definition 3.9: (Liu 2010, 2011) Let $\left\{C_t,t\ge 0\right\}$be an uncertain process.

1. $C_0=0$
2. $\left\{C_t,t\ge 0\right\}$

3. every increment $C_{t+s}-C_s$is a normal uncertainty variableb with expected value 0 and variance$t^2$, i.e., the uncertainty distribution of $C_{t+s}-C_s$is

$${\Psi }_{C_{t+s}-C_s}\left(z\right)={\left(1+\exp \left(-\frac{xz}{\sqrt{3}t}\right)\right)}^{-1}$$E12

Then $\left\{C_t,t\ge 0\right\}$is called a canonical process.

Remark 3.10: Comparing to Brownian motion process $\left\{B_t,t\ge 0\right\}$in probability theory, which is continuous almost everywhere and nowhere is differentiable, while Liu's canonical process $\left\{C_t,t\ge 0\right\}$is Lipschitz-continuous and if$\left\{C_t,t\ge 0\right\}$is differentiable somewhere, the derivative is bounded. Therefore$\left\{C_t,t\ge 0\right\}$is smoother than$\left\{B_t,t\ge 0\right\}$.

Definition 3.11: (Liu, 2010, 2011) Suppose $\left\{C_t,t\ge 0\right\}$is a canonical process, and $f$and $g$are some given functions, then

$$d{\xi }_t=f\left(t,{\xi }_t\right)dt+g\left(t,{\xi }_t\right)d{C}_t$$E13

is called an uncertain differential equation. A solution to the uncertain differential equation is the uncertain process$\left\{{\xi }_t,t\ge 0\right\}$satisfying it for any $t>0.$

Remark 3.12: Since $d{C}_t$ and $d{\xi }_t$ are only meaningful under the umbrella of uncertain integral, i.e., the an uncertain differential equation is an alternative representation of

$${\xi }_t={\xi }_0+{\displaystyle \underset{0}{\overset{t}{\int }}f\left(s,{\xi }_s\right)ds}+{\displaystyle \underset{0}{\overset{t}{\int }}g\left(s,{\xi }_s\right)d{C}_s}$$E14

Definition 3.13: The geometric canonical process $\left\{G_t,t\ge 0\right\}$satisfies the uncertain differential equation

$$d{G}_t=\alpha G_{t}dt+\sigma G_{t}d{C}_t$$E15

has a solution

$$G_t=\exp \left(\alpha t+\sigma C_t\right)$$E16

where $\alpha$can be called the drift coefficient and $\sigma >0$can be called the diffusion coefficient of the geometric canonical process $\left\{G_t,t\ge 0\right\}$due to the roles played respectively.

4. Spatial interpolation and extrapolation via inverse distance approach

Statistically, spatial interpolation and extrapolation modeling is actually a kind of linear regression modeling exercises, say, kriging methodology. Considering the shortage of California PM₁₀ data records, we will utilize a weighted linear combination approach, which was first proposed by Shepard (1968). The weights are the inverse distances between the missing value cell to the actual observed PM₁₀ value cells. The weight construction is a deterministic method, which is neutral and does not link to any specific measure theory. It is widely used in spatial predictions and map constructions in geostatistics, but is not probability oriented, rather, molecular mechanics stimulated. A unique aspect of geostatistics is the use of regionalized variables which are variables that fall between random variables and completely deterministic variables. The weight of an observed PM₁₀ value is inversely proportional to its distance from the estimated value.

$$C_{ij}=\frac{w_{ij}}{{\displaystyle \sum _{i=1}^{213}{w}_{ij}}}{z}_{ij}$$E17

Figure 4.

The 7 sites from the selected 7 counties with completed 19 year observations of PM₁₀

We wrote a VBA Macro to facilitate the interpolations and the extrapolations to "fill" up the 2048 missing value cells in terms of the 1639 cells with PM₁₀ values. With the interpolations and the extrapolations, every site has 19 PM₁₀ values now. As to whether the inverse distance approach can facilitate highly accurate predictions for each cell without a observed PM₁₀ value, we performed a re-interpolation and re- extrapolation scheme (by deleting a true PM₁₀ record, then fill it by the remaining records one by one) to evaluate the mean square value for error evaluation, the calculated mean of sum of error squares is 59.885, which is statistically significant (asymptotically).

We plotted sites 2045, 2744, 2199, 2263, 2297, 2914, and 2248 (appeared in Fig. 3) respectively in Fig. 4. By comparing Fig. 3 and 4, it is obvious that only Site 2744 the hazard level changed (moving up to next higher hazard level), while the hazard level of other six sites are unchanged. This may give an justification of the inverse distance approach. Keep in mind, the aim of this article is investigate whether the PM₁₀ level is changed over 1989 to 2007 19-year period. The change is not necessarily be accurate but reasonably calculated because of the impreciseness features of PM₁₀ complete records.

5. Uncertain analysis of site temporal pattern

Once the interpolations and the extrapolations in terms of the inverse distance approach is completed, a "complete" data set is available, containing 4047 data records of 213 sites over 19 years. The next task is for a given site, how to model the uncertain temporal pattern. It is obvious that the "complete" data set contains impreciseness uncertainty due to the interpolations and the extrapolations. We are unsure that the impreciseness uncertainty is of random uncertainty, so that we still use uncertain measure theory to pursue the temporal uncertainty modelling.

Recall that the Definition 3.13 in Section 3 facilitates a uncertain geometric canonical process,$\left\{G_t,t\ge 0\right\}$. Notice that $G_0=0$may not fit the data reality so that we propose a modified uncertain geometric canonical process, $\left\{G_t^{*},t\ge 0\right\}$with$G_0>0$:

$$G_t^{*}=G_0{G}_t=G_0\exp \left(\alpha t+\sigma C_t\right)$$E18

Note that

$$\ln G_t^{*}=\ln G_0+\alpha t+\sigma C_t$$E19

Let$y_t=\ln G_t^{*}$, ${\alpha }_0=\ln G_0$, then we have

$$y_t={\alpha }_0+\alpha t+\sigma C_t$$E20

Recall the relevant definitions in Section 3, we have

$${\rm E}\left[C_t\right]=0,\text{ and} V\left[C_t\right]=t^2$$E21

But note that for$\forall s<t$,

$$\begin{array}{l}{\rm E}\left[C_t{C}_s\right]={\rm E}\left[\left(C_s+\left(C_t-C_s\right)\right)C_s\right]\\ ={\rm E}\left[C_s^2\right]+{\rm E}\left[C_s\left(C_t-C_s\right)\right]\\ =s^2+{\rm E}\left[C_s\left(C_t-C_s\right)\right]\end{array}$$E22

Notice that the increment$C_t-C_s$is independent of$C_s$, i.e., $C_{t-s}$is independent of$C_s$.

Therefore,

$$\begin{array}{l}{\rm E}\left[C_{t-s}{C}_s\right]={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{z}_1{z}_{2}d{\Phi }_{C_{t-s},C_s}\left(z_1,z_2\right)}}\\ ={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{z}_1{z}_{2}d\min \left\{{\Psi }_{C_{t-s}}\left(z_1\right),{\Upsilon }_{C_s}\left(z_2\right)\right\}}}\\ ={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{z}_1{z}_{2}d\min \left\{{\left(1+\exp \left(-\frac{\pi z_1^2}{\sqrt{3}\left(t-s\right)}\right)\right)}^{-1},{\left(1+\exp \left(-\frac{\pi z_2^2}{\sqrt{3}s}\right)\right)}^{-1}\right\}}}\end{array}$$E23

since if ${\xi }_1$and ${\xi }_2$are independent uncertain variables with uncertainty distributions ${\Psi }_{{\xi }_1}$and ${\Upsilon }_{{\xi }_2}$ respectively, then the joint uncertainty distribution of $\left({\xi }_1,{\xi }_2\right)$ is${\Phi }_{{\xi }_1,{\xi }_2}\left(z_1,z_2\right)=\min \left\{{\Psi }_{{\xi }_1}\left(z_1\right),{\Upsilon }_{{\xi }_2}\left(z_2\right)\right\}$. Hence we obtain the expression of${\sigma }_{s,t}$:

$${\sigma }_{s,t}={\rm E}\left[C_s{C}_t\right]=s^2+{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{z}_1{z}_{2}d\min \left\{{\left(1+\exp \left(-\frac{\pi z_1^2}{\sqrt{3}\left(t-s\right)}\right)\right)}^{-1},{\left(1+\exp \left(-\frac{\pi z_2^2}{\sqrt{3}s}\right)\right)}^{-1}\right\}}}$$E24

Then the $i^{th}$"variance-covariance" matrix for uncertain vector$\left(y_{i1},y_{i2},\cdots ,y_{i19}\right)\text{'}$

$${\Gamma }_i={\sigma }^2{\left({\sigma }_{j,k}^i\right)}_{19\times 19}$$E25

where $i$is the site index, $j,k=1,2,\cdots ,19$are the entry pair in${\Gamma }_i$matrix. Hence we have a regression model (Draper and Smith, 1981, Guo et al., 2010, Guo, 2010).

For the $i^{th}$site, the regression model is

$$y_{it}={\alpha }_{i0}+{\alpha }_{i}t+{\sigma }_i{C}_{i,t}$$E26

Then in terms of the weighted least square criterion we can define an objective function as

$$J_i({\alpha }_{i0},{\alpha }_i)={\left({\underset{¯}{Y}}_i-X{\underset{¯}{\beta }}_i\right)}^{\text{'}}{\Gamma }^{-1}\left({\underset{¯}{Y}}_i-X{\underset{¯}{\beta }}_i\right)$$E27

where

$${\underset{¯}{Y}}_i=\left[\begin{array}{c}{y}_{i1}\\ y_{i2}\\ ⋮\\ y_{i19}\end{array}\right]\underset{¯}{\beta }=\left[\begin{array}{c}{\alpha }_{i0}\\ {\alpha }_i\end{array}\right]X=\left[\begin{array}{cc}\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}& \begin{array}{c}1\\ 2\\ ⋮\\ 19\end{array}\end{array}\right]$$E28

We further notice that

$$\begin{array}{l}{r}_{ij}=y_{i,j}-y_{i,j-1}\\ ={\alpha }_i+{\sigma }_i\left(C_{ij}-C_{i,j-1}\right)={\alpha }_i+{\sigma }_i\Delta C_{ij}\\ j=2,3,\cdots ,19\end{array}$$E29

Then it is reasonable to estimate ${\alpha }_i$ by

$${\widehat{\alpha }}_i=\frac{1}{18}{\displaystyle \sum _{j=2}^{19}{r}_{ij}}$$E30

Furthermore, we notice that

$${\widehat{\sigma }}_i=\sqrt{\frac{1}{18}{\displaystyle \sum _{j=2}^{19}{\left(r_{ij}-{\widehat{\alpha }}_i\right)}^2}}$$E31

Also, we can evaluate

$${\sigma }_{j,k}^i=j^2+{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{z}_1{z}_{2}d\min \left\{{\left(1+\exp \left(-\frac{\pi z_1}{\sqrt{3}\left(k-j\right)}\right)\right)}^{-1},{\left(1+\exp \left(-\frac{\pi z_2}{\sqrt{3}j}\right)\right)}^{-1}\right\}}}$$E32

in terms of numerical integration, Then an estimate for ${\Gamma }_i$matrix is obtained:

$${\widehat{\Gamma }}_i={\widehat{\sigma }}_i^2{\left({\sigma }_{j,k}^i\right)}_{19\times 19}$$E33

Finally, we use the approximated objective function

$${\widehat{J}}_i({\alpha }_{i0},{\alpha }_i)={\left({\underset{¯}{Y}}_i-X{\underset{¯}{\beta }}_i\right)}^{\text{'}}{\widehat{\Gamma }}_i^{-1}\left({\underset{¯}{Y}}_i-X{\underset{¯}{\beta }}_i\right)$$E34

to obtain a pair of estimates$\left({\tilde{\alpha }}_{i0},{\tilde{\alpha }}_i\right)$. Repeat this estimation process until all the 213 weighted least square estimate $\left({\tilde{\alpha }}_{i0},{\tilde{\alpha }}_i\right)$are obtained.

Recall the definition of coefficient ${\alpha }_i$ so that the sign and the absolute value of ${\alpha }_i$ indicates the geometric change over the 19 years. Since the estimation procedure of ${\alpha }_i$ involves all the spatial-temporal information, it is reasonable to have them plotted in a kriging map to reveal the overall changes over 19-year period.

6. Kriging maps and time-change maps based on completed PM₁₀ data

Kriging map presentation is vital for a geostatistian's visualization, and maps reveal hidden information or the whole picture. A sample statistic is typically condensing the wide-spread information into a numerical point. While, a kringing map is actually a map statistic (or a statistical map) which contains infinitely many information aggregated from limited "sample" information (i.e., observations). Kriging itself is not specifically probability oriented, it is another weighted linear combination prediction, but requires more mathematical assumptions. In fuzzy geostatistics, the fuzzy kriging scheme has also been developed (Bardossy et al., 1990).

Ordinary kriging (abbreviated as OK) is a linear predictor, see Cressie (1991) and Mase (2011). The formula is

$$Z\left(s_0\right)={\displaystyle \sum _{j=1}^N{\lambda }_{j}Z\left(s_j\right)}$$E35

where $s_j$ are spatial locations with observation $Z\left(s_j\right)$available and the coefficients ${\lambda }_j$ satisfy the OK linear equation system

$$\{\begin{array}{c}{\displaystyle \sum _{j=1}^N{\lambda }_j\gamma \left(\epsilon \left(s_j\right)-\epsilon \left(s_i\right)\right)}-\psi =\gamma \left(\epsilon \left(s_0\right)-\epsilon \left(s_i\right)\right),i=1,2,\cdots ,N\\ {\displaystyle \sum _{j=1}^N{\lambda }_j}=1\end{array}$$E36

The OK system is generated under the assumptions of an additive spatial model

$$Z\left(s\right)=\mu \left(s\right)+\epsilon \left(s\right)$$E37

where µ(s) is the basic (expected) spatial trend and $\epsilon \left(s\right)$ is a Gaussian error$N\left(0,{\sigma }^2\left(s\right)\right)$, i.e., Gaussian variable with mean and variance

$${\rm E}\left[\epsilon \left(s\right)\right]=0,V\left[\epsilon \left(s\right)\right]={\sigma }^2\left(s\right)$$E38

respectively. Accordingly, the variogram $2\gamma$ of the random error function $\epsilon (\cdot )$ is just defined by

$$2\gamma \left(h\right)={\rm E}\left[{\left(\epsilon \left(s+h\right)-\epsilon \left(h\right)\right)}^2\right]$$E39

where $h$ is the separate vector between two spatial point $s+h$ and $s$under the isotropy assumption.

Figure 5.

Kriging Prediction Maps for PM₁₀ in California 1989-2007.

The 213 observation sites now have 19-year PM₁₀ values, a "complete" data set is now available, containing 4047 data records of 213 sites over 19 years, and then the 19 ordinary kriging pred4iction maps are generated for comparisons. In Fig. 5, all 19 years of PM₁₀ concentration in California State are shown. It is very interesting to examine the change in PM₁₀ concentrations through the 19 years, based upon the modelled complete 213 site data. In particular, 1998 shows to have an extremely low PM₁₀ concentration. Although air quality is varied over the years, but in general, the PM₁₀ concentration is decreasing, showing an improvement of air quality trend.

Figure 6.

Changes in PM₁₀ values and the rate of change of PM₁₀ in California between 1989 and 2007.

As one can clearly see from Fig. 6, that PM₁₀ concentration has clearly decreased over the 19 years, and air quality has improved remarkably over the years. The blue and green colours show negative changes, and red shows positive changes or near positive changes. Counties such as San Diego, Inyo, Santa Barbara, Imperial, still show an increase in PM₁₀ concentration in the air, and indicate bad air quality. While Kern, Modoc, Siskiyou counties show the most improvement in air quality. The left map in Fig. 6 is PM₁₀ record difference between 2007 and 1989 at each location, in total 231 values, and then a difference map is constructed. It is obvious that the difference map only utilizes 1989 and 2007 two-year PM₁₀ records, 1990, 1991,..., 2006 seventeen years' information do not participate the change map construction. The right map in Fig. 6 show completed${\tilde{\alpha }}_i,i=1,2,\cdots ,213$, the rate of change over 1989 to 2007 19-year period.

Note that the calculations of ${\tilde{\alpha }}_i,i=1,2,\cdots ,213$involve all nineteen years by temporal regression, the dependent variable $y$are estimated form the actual PM₁₀ observations cross over all the available locations. Therefore, the rate of change parameter ${\alpha }_i$at each individual location contains all spatial-temporal information. It is reasonable to say the rate of change parameter ${\tilde{\alpha }}_i$is an aggregate statistic for revealing the 19-year changes over 213 locations. ${\tilde{\alpha }}_i$kriging map is thus different from 2007-1989 kriging maps. The positive sign of ${\tilde{\alpha }}_i$ indicates the increasing trend in PM₁₀ concentration, while the negative sign f ${\tilde{\alpha }}_i$ indicates the decreasing trend in PM₁₀ concentration. The absolute value of ${\tilde{\alpha }}_i$reveals the magnitude of change of PM₁₀ concentration. It is worth to report, among 213 locations, 193 locations have negative${\tilde{\alpha }}_i$, while the negative ${\tilde{\alpha }}_i$locations are 20 (9% approximately).

7. Discussion

Air quality and health is always a central issue to public concern on the quality of life. In this chapter, we examined PM₁₀ levels over 19 years, from 1989 to 2007, in the California State. Facing the difficult task of a lack of "complete" PM₁₀ observational data, we utilised the inverse distance weight methodology to "fill" in the locations with missing values. By doing so, the impreciseness uncertainty is introduced, which is not necessarily explained by probability measure foundation. We noted the character of a regionalized variable in geostatistics and therefore engage Liu's (2010, 2011) uncertainty theory to address the impreciseness uncertainty. In this case, we developed a series of uncertain measure theory founded spatial-temporal methodology, including the inverse distance scheme, the kriging scheme, and the geometric canonical process based weighted regression analysis in order to extract the change information from the incomplete 1989-2007 PM₁₀ records. The use of the rate of change parameter alpha is a new idea and it is an aggregate change index utilized all spatial-temporal data information available. It is far better than classical change treatments. However, due to the limitations of our ability, we are unable to demonstrate the detailed uncertain measure based spatial analysis model. In the future research, we plan to develop a more solid uncertain spatial prediction methodology.

Acknowledgments

I would like to thank the California Air Resources Board for providing the air quality data used in this paper. This study is supported financially by the National Research Foundation of South Africa (Ref. No. IFR2009090800013) and (Ref. No. IFR2011040400096).