Semi-Supervised Learning for Integration of Aerosol Predictions from Multiple Satellite Instruments

1 Introduction

Aerosols are small airborne particles produced by natural and man-made sources that both reflect and absorb incoming Solar radiation. Depending on their distribution and composition, aerosols can result either in cooling or warming of the atmosphere, thus having a major role in regulating the climate system. Distribution of aerosols is measured by Aerosol Optical Depth (AOD or τ), a quantitative measure of the extinction of Solar radiation by scattering and absorption between the top of the atmosphere and the surface. AOD is an important input to climate models, and it can significantly impact predictions of future climate changes [Randall et al., 2007]. Considering that climate predictions influence decisions of policy makers, accurate AOD estimation is a task of global significance. In addition to its impact on climate studies, AOD is an important quantity in estimation of air pollution. For example, it was shown in [Liu et al., 2009] that AOD is an accurate predictor of PM2.5, the concentration of particulate matter with aerodynamic diameters ≤ 2.5μm, which poses a serious health hazard to the population [Hu and Rao, 2009].

Currently, a number of instruments aboard several Earthobserving satellites report their AOD estimates, such as MODIS instrument aboard Terra and Aqua satellites [King et al., 2003], MISR aboard Terra [Diner et al., 1998], OMI aboard Aura [Torres et al., 2002], SeaWiFS aboard SeaStar [Wang et al., 2000], and others. All these instruments have a capability of providing global estimates of AOD distribution with a fine spatial (few kilometers) and temporal (few days) resolution. Each instrument has different properties and estimates AOD using a different algorithm developed by domain scientists. Coverage and quality of satellite measurements can differ from instrument to instrument for a number of reasons. As illustrated in Figure 1, width of the field of view of MODIS instrument is 2,330km, allowing MODIS to observe the entire Earth every day, as opposed to 360km width of MISR instrument, which requires 9 days for global coverage. The quality of AOD estimates from different instruments varies with atmospheric and surface conditions [Mishchenko et al., 2010]. For example, 9 cameras observing Earth at 9 different angles used by MISR allow it to be more accurate than MODIS when clouds are present, over bright surfaces, or for some types of aerosol compositions. In addition to satelliteborne sensors, AOD is also measured by a network of groundbased sensors from AERONET [Holben et al., 1998], placed at several hundred unevenly distributed locations across the globe, see Figure 2. AERONET AOD measurements are considered a ground-truth, as they are several times more accurate than the best available satellite AOD estimations. The drawback of AERONET is that it has a very limited spatial coverage, and that it cannot be used to provide global estimation of AOD distribution required for climate models.

Different spatial and temporal coverage, design, and specific mission objectives of the satellite-borne instruments mean that they observe and measure different, possibly complementary aspects of the same phenomenon. Instead of considering AOD estimates of individual instruments in isolation, combining measurements from different sources into an aggregated estimate may prove to be the best path towards obtaining a higher-quality global AOD estimation. A recent study by [Mishchenko et al., 2010] confirmed this hypothesis by illustrating that simple average of collocated Terra MODIS and MISR AOD estimations resulted in improved accuracy.

The combination of experts that ultimately yields an estimate that is more accurate than any of the individual forecasts is a well-researched topic. Assuming the Gaussian distribution of prediction errors of AOD and no missing experts, [Bates and Granger, 1969; Granger and Ramanathan, 1984] proposed how to learn the optimal combination of experts from labeled data. If data set is unlabeled, [Ristovski et al., 2010] proposed how to learn a combination of experts by extending the classification-based method from [Raykar et al., 2009]. However, the approach assumed that experts are independent, and that all experts are available for aggregation.

We propose a novel method suitable for finding a linear combination of AOD estimations from multiple instruments. There are several interesting challenges specific to the aerosol domain that had to be addressed. (1) As quality of different instruments varies with atmospheric and surface conditions, it is not likely that the same linear combination would work equally well at different locations, for example, in North America and Africa [Levy et al., 2007]. Therefore, it might be needed to develop specialized combinations for different regions around the globe. (2) Number of labeled data points is relatively small. For example, in North America, thanks to a relative abundance of AERONET sites, the number of labeled data points can exceed a thousand every year, while in Africa and parts of Asia there are very few AERONET sites, and the number of labeled data points could be measured in tens every year. In addition to their small number, labeled data points might cover only a limited set of conditions observable at AERONET locations. On the other hand, the number of unlabeled data points is orders of magnitudes larger. An open question in AOD estimation is how to exploit labeled and unlabeled data. (3) As shown in Figure 1, which illustrates daily coverage of different sensors over the USA, for most of the labeled (e.g., points A and B) and unlabeled (e.g., points C and D) data points, AOD estimations from some of the instruments are missing. For example, points A and C have AOD estimate from all 3 satellite instruments, while points B and D are just outside of MISR's field of view and do not have its AOD estimate. This opens a question of learning from data with significant amount of missing AOD estimations.

In this paper, we assume that estimation errors of individual satellite instruments have multivariate Gaussian distribution, and propose a semi-supervised method that can handle missing data while being able to partition the data into homogeneous subsets on which specialized aggregators are learned. We note that our method can be seen as a significant generalization of the traditional supervised method for combination of experts by [Bates and Granger, 1969; Granger and Ramanathan, 1984], as well as of recently proposed unsupervised method for averaging of experts in regression by [Ristovski et al., 2010].

2 Methodology

2.1 Problem setup and assumptions

Let us assume that we have a training data set D = {{ˆ yik}k=1,...,K, yi}i=1,...,N, where target value yi for the ith data point is predicted by K experts, with the kth expert providing an opinion in a form of prediction ˆ yik. For example, in the aerosol domain that we study, the experts are satellite instruments and the predictions are their individual AOD estimates. We assume that data points are independent and identically distributed (IID), and that ground truth yi is normally distributed with mean μy and variance σ2 y ,

We also assume that the first Nu data points are unlabeled, while the last Nl data points are labeled, i.e., we have a ground truth only for data points indexed by i = (Nu + 1), . . . ,N, with N = (Nu + Nl). We use 1 to denote a column-vector of all ones, 0 to denote a matrix of all zeros, and ˆyi = [ˆ yi1, . . . , ˆ yiK]T to denote a column-vector of expert predictions for the ith data point. We assume expert predictions for the ith data point are sampled from a multivariate normal distribution as

This assumption allows the experts to be correlated (i.e., Σ is non-diagonal), as is the case in aerosol domain. We first consider a case where all experts are available, and then extend the methodology to account for missing experts. Given D, the objective is to learn Σ, μy, and σ2 y . By Θ = {Σ, μy, σ2 y } we denote a set of parameters to be learned.

Once Θ is learned, and given expert predictions ˆyi, aggregated prediction yi for the ith data point can be found as a mean of the posterior distribution yi|ˆyi ∼ N(yi, (1TΣ′-11)-1 ), where mean yi is computed as

ˆy′i = [ˆyT i , μy]T, and Σ′ is a (K+1)×(K+1) block matrix

2.2 Semi-supervised combination of experts

Given the parameters Θ of the model, the probability of observing the data set D can be written as

where subscripts u and l denote unlabeled and labeled parts of the data set, respectively. Let us first consider P(Du|Θ). As the data points are sampled IID, the probability factorizes over individual data points, and we can write

As both probabilities under the integral are assumed Gaussian, by solving the integral we obtain

Prior parameters μy and σ2 y , appearing in ˆy′i and Σ′ , respectively, can be fixed to some values, such as mean and variance of the available target values, or could be learned. In order to keep the notation simple, in the remainder of the section we assume σ2 y → ∞, which amounts to an uninformative prior P(y|Θ). We note that it is straightforward to modify the following expressions for finite σ2 y , or to derive a learning rule.

Likelihood of the labeled part can be written as

which, due to (2), is a product of Nl multivariate Gaussians. Then, combining equations (5), (7), and (8), we can compute the likelihood of the data set D. After finding the derivative of the log-likelihood with respect to Σ-1 and equating the resulting expression with zero, we obtain the following expression for computing Σ matrix,

where ˆYu and ˆYl are Nu×K and Nl ×K matrices of expert predictions for unlabeled and labeled data, respectively, with each row corresponding to a single data point, and yl is an Nl × 1 column-vector of ground-truth values. Equation (9) yields an iterative procedure for learning Σ, where Σ on the l.h.s. is a new value, and Σ on the r.h.s. is an old value.

2.3 Missing experts

Let us now consider the case where some experts are missing. For example, let us assume that the ith data point has q missing predictions. Then, we reorganize vector ˆyi in such a way so that the first a = (K - q) elements are available predictions, while the last q elements are missing predictions, i.e., ˆyi = [ˆyT ai , ˆyT qi ]T. Similarly, we reorganize Σ-1 matrix so that the first a rows/columns correspond to available predictions, while the remaining q rows/columns correspond to missing predictions, or

where Πi is a permutation function used to reorder both rows and columns of Σ-1 according to the ith data point, and U is an a × a matrix. Given the covariance matrix Σ and a vector of expert predictions ˆyai for the ith data point, the aggregated prediction yi can be found as a mean of the posterior distribution yi|ˆyai ∼ N(yi, (1TU′i1)-1 ), where we introduced U′ = U-VQ-1VT to simplify the notation, and

Note that we appended subscript i to indicate that the size of a matrix U′i depends on the number of available experts for the ith data point.

In the following, we derive the update equation for Σ. The probability of observing the ith unlabeled point is equal to

Solving the equation (12) we obtain

In a very similar manner we can find the probability of observing the ith labeled data point. It follows

which, after solving the integral, results in

By combining equations (5), (13), and (15), we can find the likelihood of the data set D. After finding derivative of the log-likelihood with respect to Σ-1 [Brewer, 1978] and equating the resulting expression with zero, we obtain the following expression for computing Σ matrix

where Π-1 i is an inverse permutation function that reorders rows and columns of the matrix back to the original order of experts, symmetric K ×K matrix Ψi is equal to

and llbracketAirrbracket for some symmetric a × a matrix Ai denotes the following symmetric K × K matrix

2.4 Incorporating prior probability P(Θ)

Let us consider the case where we have some prior knowledge about Σ, and would like to include this knowledge into the model. As we assumed σ2 y → ∞, it follows Θ = {Σ}, and we can write

Note that we defined prior P(Σ-1 ) in terms of an inverse of the covariance matrix (i.e., in terms of a precision matrix). For the precision matrix Σ-1 we choose the prior as a Wishart distribution W(S, n) with given K × K scale matrix S and n > (K - 1) degrees of freedom, resulting in

which is a conjugate prior for multivariate Gaussian distribution, and where ΓK is the multivariate gamma function. After choosing n = (K + 2) and finding the derivative of the log-likelihood with respect to Σ-1 , we obtain the following update equation for the covariance matrix Σ

2.5 Data partitioning using a latent variable

It is an inherent property of the experts in the aerosol domain that they do not maintain the same quality of predictions across all observed conditions. To address this characteristic of the aggregation problem, we consider partitioning the data points into several groups, called the regimes, where each regime is governed by a different multivariate Gaussian from (2) [Weigend et al., 1995]. In the following we assume there are R regimes, and that we have available a feature vector xi for the ith data point that could be used to assign it to an appropriate regime.

Assuming a mixture of R regimes, probability of observing expert predictions ˆyai for the ith labeled data point can be written as

where Pr(ˆyai|yi) = P(ˆyai|regimer , xi, yi,Θ), πir(xi) = P(regimer |xi,Θ), and where appended subscript r denotes the rth regime. Similarly, we can write probability of observing expert predictions ˆyai for the ith unlabeled data point as

Probability of observing the ith unlabeled or labeled data point given that it was generated by the rth regime, Pr(ˆyai) or Pr(ˆyai|yi), respectively, can be computed by considering equations (13) and (15), respectively. The aggregated prediction yi can be found as

To facilitate model optimization, we consider regime assignments as unobserved data, and introduce a latent indicator variable zir such that

Further, by introducing zi = [zi1, . . . , ziR]T, we can write the complete-data likelihood for the ith labeled data point as

Note that, due to the lack of space, in equations (26), (27), and (29), we only give expressions for labeled data. However, when dealing with the ith unlabeled data point we simply need to replace Pr(ˆyai|yi) by Pr(ˆyai). Then, the complete-data log-likelihood L is equal to

Expectation-Maximization (EM) algorithm [Dempster et al., 1977] can be used to find the parameters Θ that maximize L from (27).

EM algorithm for semi-supervised aggregation

Before moving on, we need to decide on the parameterization of the prior probability πir. We define this probability using a softmax function,

where we defined a prototype vector qr and feature scaling matrix Λr for each regime, to be found during optimization, resulting in Θ = {Σr, qr, Λr}r=1,...,R.

In the E-step, we compute the current expectation of posterior probability hir that the rth regime is "responsible" for generating expert predictions for the ith labeled data point as

Then, in the M-step, we fix values of hir for all data points and regimes, and optimize L with respect to covariance matrices Σr and prototype vectors and scaling matrices qr, Λr, r = 1, . . . , R. Note that the derivatives of L with respect to these two sets of variables are independent from each other, and the optimization of Σr on one side, and qr and Λr on the other, can be easily parallelized. After derivation, the update equation for Σr can be written as

while prototype vector qr and scaling matrix Λr are found through the gradient ascent optimization as

where η is an appropriately set learning rate.

3 Experiments

In this section, we first experimentally validate the semisupervised aggregation on synthetic data, and then apply the method to AOD estimation using real-world aerosol data set.

3.1 Validation on synthetic data

We started by evaluating our method on synthetic data generated as follows: for a given number of regimes R and experts K, we selected a prototype and a covariance matrix for each regime. Then, we assigned the ith data point uniformly at random with probability 1/R to a regime, say the lth regime, and obtained features xi by sampling frommultivariate Gaussian with mean ql and diagonal covariance matrix 0.5I. We sampled ground-truth value yi from zero-mean Gaussian with unit-variance, then sampled K expert predictions from a Gaussian N(yi1, Σl). Finally, we removed each expert's prediction with probability 0.5 to simulate missing experts. In all experiments we set S = I, and used 15 EM iterations. Learning rate η was set through cross-validation.

First, in order to evaluate the semi-supervised method without clustering, we set K = 5, R = 1, and Σ1 = diag([0.1, 0.2, 0.3, 0.4, 0.5]). We compared our learning method to a baseline method that averages all available experts, as well as to the optimal predictor that computes the prediction (24) using the true Σ1. We increased the number of training points N from 10 to 100 in increments of 10, and for each N we experimented with percentage of labeled data points equal to 0%, 20%, 40%, and 100% (shown as four solid lines in Figure 3). The results in terms of Root Mean Squared Error (RMSE), evaluated on 1,000 testing points generated in the same way as the training set and averaged over 100 experiments, are shown in Figure 3a. We can see that the performance of the fully unsupervised approach, given by the top-most full line, is already better than simple averaging, which further improves as the number of unlabeled data grows. Moreover, as we increase the number of labeled points, the semi-supervised method further improves the accuracy, approaching the lower bound on RMSE achieved by the optimal combination of experts.

Next, we generated the data using two regimes by setting q1 = [1, 1], q2 = [-1,-1], Σ1 = diag([0.1, 0.2, 0.3, 0.4, 0.5]), Σ2 = diag([0.5, 0.4, 0.3, 0.2, 0.1]), and we set R = 2. The results in terms of RMSE are given in Figure 3b, where we also show accuracy of the proposed method which used only labeled data, but assumed only one cluster. The RMSE of supervised method that assumed only a single cluster is worse than simple averaging, and approached it as the data size increased. Unsupervised method using two clusters achieved better accuracy than simple averaging, and RMSE further decreased with larger data sizes. Introduction of labeled data points further decreased the RMSE.

3.2 Validation on aerosol data

To construct an aerosol data set we considered ground-based AERONET data 1 from the United States and collocated data from 5 satellite instruments 2 spanning years 2006 to 2010. We considered AERONET measurements at 10:30am local time as the ground truth. Among the K = 5 instruments, two measured AOD at around 10:30am local time (Terra MODIS, MISR), and three at around 1:30pm local time (Aqua MODIS, OMI, SeaWiFS). After removing AERONET sites with too few observations, there remained 33 sites in the data set, with locations shown in Figure 4. This resulted in a labeled data set with N = 6,913 data points, where 58% of expert predictions were missing. We used this data set for two sets of experiments: (1) evaluating usefulness of partitioning; and (2) evaluating usefulness of unlabeled data.

In both sets of experiments we performed leave-one-siteout cross-validation. In the first set of experiments, from each training site we randomly sampled 100 data points, and assumed that 50 of them are labeled and 50 unlabeled. From the left-out site we sampled 100 data points. We used the geographic location (i.e., longitude and latitude) of the corresponding AERONET site as a feature vector xi for the ith data point. We used our proposed method with R = 1 and R = 2, repeating the experiments 5 times. We compared the performance to a baseline method that takes a simple average of available expert predictions. RMSE is reported at the top of Table 1. We can see that semi-supervised aggregation for both R = 1 and R = 2 had significantly lower RMSE than the baseline. Moreover, by increasing the number of clusters from 1 (i.e., without clustering) to 2 we observed a drop in RMSE of nearly 5%. In Figure 4 we color-code the AERONET sites according to their cluster assignments. Interestingly, the clustering roughly corresponds to partitioning proposed by domain scientists [Levy et al., 2007]. For the south-western cluster the weights of linear combination assigned to MISR, Terra MODIS, Aqua MODIS, OMI, and SeaWiFS instruments, given predictions of all the experts, were [0.51, 0.31, 0.09, 0.01, 0.08], while for the north-eastern cluster they were [0.24, 0.27, 0.21, 0.15, 0.12], respectively. Consistent with the domain knowledge, MISR obtained the largest weight in the first cluster, while all instruments were given similar weights in the second cluster.

1 aeronet.gsfc.nasa.gov/cgi-bin/combined data access new

2 disc.sci.gsfc.nasa.gov/aerosols/services/mapss/mapssdoc

In the second set of experiments, we simulated conditions consistent with aerosol data availability in Africa and large parts of Asia, where very few AERONET sites are available and labeled data are very scarce. We randomly selected 2, 4, or 6 training AERONET sites and took 100 data points from each of them as labeled data. Then, we selected 100 data points from the remaining training AERONET sites and treated them as unlabeled data. We trained one model which used only labeled data, and one that used both labeled and unlabeled data. We used R = 2 clusters in both cases. The results given at the bottom of Table 1 show that the RMSE of purely supervised approach decreased with the number of AERONET sites. More importantly, unlabeled data were helpful and led to significant reductions in RMSE. This benefit of unlabeled data increased with the decrease in the number of labeled data points. The results in Table 1 confirm the validity of the proposed semi-supervised method for aggregation of experts, which is able to account for missing experts, find a partition of data into clusters, and construct specialized aggregators on each cluster.

4 Conclusion

We proposed a semi-supervised method for aggregation of AOD predictions from noisy satellite-borne sensors into a single, more accurate estimate. By assuming that expert predictions follow multivariate Gaussian distribution, the method accounts for both missing experts and unlabeled data in a principled manner, addressing an issue inherent to the remote sensing domain. Moreover, we also cluster the data during training by introducing a latent indicator variable for each cluster, resulting in a more interpretable model. Results on synthetic and real-world aerosol data comprising 5 satelliteborne sensors indicate the benefits of the proposed approach.

Acknowledgments

This work was supported by NSF grant IIS-1117433.