Calibrated sky imager for aerosol optical properties determination

Calibrated sky imager for aerosol optical properties determination A. Cazorla, J. E. Shields, M. E. Karr, A. Burden, F. J. Olmo, and L. Alados-Arboledas Departamento de Fı́sica Aplicada, Facultad de Ciencias, Universidad de Granada, Fuentenueva s/n, 18071, Granada, Spain Centro Andaluz de Medio Ambiente (CEAMA), Junta de Andalucı́a-Universidad de Granada, Avda. del Mediterraneo s/n. 18071, Granada, Spain Marine Physical Lab, Scripps Institution of Oceanography, University of California San Diego, 9500 Gilman Dr., La Jolla CA 92093-0701, USA Received: 27 August 2008 – Accepted: 7 October 2008 – Published: 28 November 2008 Correspondence to: L. Alados-Arboledas (alados@ugr.es) Published by Copernicus Publications on behalf of the European Geosciences Union.

In this sense long range transport events like Saharan dust outbreaks (Lyamani et al., 2005;Lyamani et al., 2006a, b) or global scale events like stratospheric aerosols following major volcanic eruptions like El Chichón and Mount Pinatubo (Olmo and Alados-Arboledas, 1995) represent extreme cases of this variability. Remote sensing appears to be a valuable tool for characterizing the physical and optical properties of the aerosol. 15 Sun-photometry is the most common way to characterize aerosol in daytime from the ground. An instrument such as a CIMEL CE318 photometer gives useful information about column-integrated physical and optical properties of the atmospheric aerosol (Holben et al., 1998). Sun-photometer networks such as Aerosol Robotic Network (AERONET) (Holben et al., 1998) address the spatial problem in aerosol characteriza-20 tion, but temporal resolution remains a problem. Another problem is the assessment of valid measurements, i.e. the CIMEL CE318 requires that the sun not be obscured by clouds, and sometimes it is difficult to perform the cloud rejection test (Smirnov et al., 2000).
Ground-based sky imagery has been used for years for cloud cover assessment e.g.
characteristics from sky imagers, we can supplement the existing aerosol data bases. In this paper, the Whole Sky Imager (WSI), a calibrated ground-based sky imager, developed in the Atmospheric Optics Group (Marine Physical Laboratory, Scripps Institution of Oceanography, University of California San Diego) (AOG) has been tested to determine optical properties of the atmospheric aerosol. Different neural network-5 based models estimate the aerosol optical depth for three wavelengths using the radiance extracted from the principal plane of sky images from the WSI (i.e. at constant azimuth angle equal to the solar azimuthal angle, with varied zenithal angles) as input parameters. TheÅngström coefficients α and β (Ångström, 1964) are also derived from the aerosol optical depth estimated with the models and a neural network-based 10 model also estimates theÅngström exponent α. The models use data from a CIMEL CE318 photometer (Holben et al., 1998) for training and validation.

Experimental site
This work is the result of cooperation between the Atmospheric Physics Group 15 (Dept. Applied Physics, University of Granada) and the Atmospheric Optics Group (Marine Physical Laboratory, Scripps Institution of Oceanography, University of California San Diego). This group has been researching and developing sky imagers for decades (Johnson et al., 1989;Shields et al., 1993Shields et al., , 1998b. They have different sky imagers in different locations. The WSI used in this work was located in the Southern general circulation models used for climate research. The SGP site consists of in situ and remote-sensing instrument clusters arrayed across approximately 143 000 square kilometers in north-central Oklahoma. Figure 1 shows a map of the facility. The central facility is a heavily instrumented location on 0.65 square kilometers of cattle pasture and wheat fields southeast of Lamont, Oklahoma (36.61 • N, 97.5 • W, 320 m a.s.l.).

5
More than 30 instrument clusters have been placed around the SGP site, at the Central Facility and at Boundary, Extended, and Intermediate Facilities. The locations for the instruments were chosen so that the measurements reflect conditions over the typical distribution of land uses within the site.
Both instruments, the CIMEL CE318 and the WSI are located in the central facility (Shields, 1998a).

The sun-photometer CIMEL CE-318
Sun photometry is the most widely used technique for atmospheric aerosol characterization in daytime. The CIMEL CE318 automatic sun tracking photometer (Holben et al., 1998) has been designed to measure sun and sky radiance in order to derive to-15 tal column water vapor, ozone and aerosol properties using a combination of spectral filters and azimuth/zenith viewing controlled by a microprocessor. The CIMEL CE318 is the standard instrument in the AERONET network (Holben et al., 1998). AERONET collaboration provides globally distributed observations of spectral AOD, inversion products (Dubovik et al., 2002(Dubovik et al., , 2006 and precipitable water in 20 diverse aerosol regimes. AOD data are computed for three data quality levels: Level 1.0 (unscreened), Level 1.5 (cloud-screened), and Level 2.0 (cloud screened and qualityassured).
The CIMEL CE318 used in this work has operated in the AERONET program since 1994. We use the level 2.0 data (quality-assured) and the parameters extracted are 2.3 The Whole Sky Imager The Atmospheric Optics Group has been very active in the development of digital sky imager for over two decades. The original concept for the WSI evolved from a measurement and modeling program using multiple sensors for monitoring sky radiance, atmospheric scattering coefficient profiles and other parameters related to vision through the atmosphere (Johnson et al., 1980(Johnson et al., , 1989. With the use of very low noise 16 bit CCD cameras and an occultor designed to handle both sun and moon, these systems were further developed into the Day/Night WSI (Shields et al., 1993(Shields et al., , 1998a(Shields et al., , 1998b(Shields et al., , and 2003. The Day/Night WSI is a 16-bit digital imaging system that acquires images of the 10 full sky (2π hemisphere) under both day and night conditions in order to assess cloud fraction, cloud morphology, and radiance distribution. The WSI measures the sky radiance in approximately 185 000 directions simultaneously by using a 512×512 CCD sensor. The result is a 34 µsteradian field of view (FOV) in each direction, to cover the full 2π steradian dome. Images are acquired in visible and near infrared (NIR) wave- 15 bands with filters at 450 nm (blue), 650 nm (red), and 800 nm (NIR) under sunlight or moonlight. Open hole is used for starlight and most moonlight conditions. The FWHM of the filters is 70 nm. A picture of the instrument fielded at the Oklahoma Cloud and Radiation Testbed (CART) site at SGP is shown in Fig. 2. The primary features seen in this figure are the environmental housing that protects the sensor and electronics, 20 the optical dome, and the solar occultor that shades the optics. Even though the camera would not be damaged by direct sun radiation, this shading is desirable because it minimizes stray light especially near the sun. Figure 3 shows a daytime image. The center of the image is the zenith and the edges are the horizon. One of the capabilities of the WSI is the determination of the absolute sky radiance 25 distribution. The fisheye lens directs the light from different directions onto different pixels in the image plane, and the signal of each pixel may be calibrated to yield a determination of the absolute sky radiance, in W/m 2 µm sr, in that direction. The FOV for 19994 each pixel is approximately 34 µsteradians. Thus, this radiance product is equivalent conceptually to a radiance distribution determined by a scanning radiometer, except that all radiances are acquired simultaneously (at 185 000 points) and at a very high spatial resolution (Shields et al., 1998b). The instrument follows a process to get the radiance and geometric calibration 5 (Shields et al., 1998b).

Methodology
The sky radiance depends directly on the aerosol load through several parameters connected with extensive and intensive aerosol properties. While previous investigations have related AOD to radiances measured in a restricted range of scattering angle 10 to simulate the spaceborne point of view (e.g. Kaufman, 1993;Sánchez et al., 1998), here the development is focused on surface measurements. There is a dependency between radiance along the principal plane and the aerosol optical parameters (Olmo et al., 2008). Radiance along the principal plane is affected by the particle size, which in turn has impacts on the α parameter and the turbidity 15 factor (β parameter). The impact of the α parameter and β parameter on the radiance are shown in Fig. 4a and b. Both parameters are directly related with the aerosol optical depth (Ångström, 1964).
These previous works, along with the results obtained with the All-Sky Imager characterizing the atmospheric aerosol (Cazorla et al., 2008b) are the basis of this work. 20 We also have considered the use of sky radiance in the principal plane since the obstruction of the image due to the shadow system is smaller.
The data set selected in this work comprises the period from 1 October 2001 to 29 September 2002. This data set is from a whole year so we can model the seasonal variability of the atmospheric aerosol. Using the cloud decision images processed by 25 AOG we sorted out all the cases with clouds, to work with the clear-sky results. A total of 1047 clear-sky image sets (i.e. 3 spectral images acquired in one set) were Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version Interactive Discussion associated with a synchronous CIMEL measurement, applying a ±5 min margin. This image set has been used to create and validate the model.

Retrieving the radiance over the principal plane of the sky images
Knowing the Sun position we can locate the principal plane making the azimuthal angle equal to the Sun azimuthal angle or that angle plus 180 degrees, and varying the zenith 5 angle. The radiance over the principal plane has been extracted for the whole data set from scattering angle 1 • to 100 • . The cloud decision image provided by the WSI has been used to apply a mask over the non valid values (horizon obstruction and shadow system obstruction). The stored values from the principal plane are calibrated values, i.e. sky radiance, in W/m 2 µm sr, for every scattering angle in the principal plane. Thus,

10
we have the sky radiance over the principal plane for scattering angles from 1 to 100 • and the Solar Zenith Angle in degrees (SZA • ) as inputs. The AOD at 440, 675 and 870 nm from the CIMEL data is used as the target data in the neural network model that we will describe in the next section.
3.2 Neural networks and radial basis function networks 15 According to Haykin (1994) a neural network resembles the human brain in two aspects: the knowledge is acquired by the network through a learning process, and interneuron connection strengths known as synaptic weights are used to store the knowledge. Once a neuron is set up, it can learn to emulate behaviors such as classification, pattern recognition, function approximation, control system, etc.

20
Neural networks have been widely used in atmospheric science recently (e.g. Gutiérrez et al., 2004) and we have experimented with them in several research applications (Alados et al., 2004(Alados et al., , 2007Cazorla et al., 2005Cazorla et al., , 2008bGil et al., 2005). The radial basis function (RBF) networks (Gutiérrez et al., 2004;Yee and Haykin, 2001) are especially suitable for function approximation. The inputs of the radial basis networks 25 are the variables of the function, and the output is the function approximation. RBF's emerged as a variant of artificial neural networks in the late 80s. However, their roots are entrenched in much older pattern recognition techniques as for example potential functions, clustering, functional approximation, spline interpolation and mixture models (Tou and Gonzalez, 1974). RBF's are embedded in a two layer neural network, where each hidden unit implements a radial activated function. The output 5 units implement a weighted sum of hidden unit outputs. The input into an RBF network is nonlinear while the output is linear. Their excellent approximation capabilities have been studied by Poggio and Girosi (1990) and Park and Sandberg (1991). Due to their nonlinear approximation properties, RBF networks are able to model complex mappings, which perceptron neural networks can only model by means of multiple intermediary layers (Haykin, 1994).
In order to use a RBF network we need to specify the hidden unit activation function, the number of processing units, a criterion for modeling a given task and a training algorithm for finding the parameters of the network. Finding the RBF weights is called network training. If we have at hand a set of input-output pairs, called a training set, we 15 optimize the network parameters in order to fit the network outputs to the given inputs. The fit is evaluated by means of a cost function, usually the mean square error. After training, the RBF network can be used with data whose underlying statistics is similar to that of the training set.
The topology of the RBF networks consists in a two-layer feed-forward neural net-20 work. Such a network is characterized by a set of inputs and a set of outputs. In between the inputs and outputs there is a layer of processing units called hidden units. Each of them implements a radial basis function. Figure 5 shows the topology of the RBF networks. RBF networks are characterized for the transfer function which is the RBF. In this 25 case input to the transfer function is the vector distance between its weight vector w and the input vector p, multiplied by the bias b The transfer function for a radial basis neuron is a Gaussian function The RBF has a maximum of 1 when its input is 0. As the distance between w and p decreases, the output increases. Thus, a radial basis neuron acts as a detector that produces 1 whenever the input p is identical to its weight vector w . The bias b allows 5 the sensitivity of the neuron to be adjusted.
This topology allows a very simple training for function approximation or interpolation. Every sample in the training set creates a new neuron in the hidden layer. The weight of the connection input-neuron is set to the input value. Therefore n in Eq. (1) is 0 and a in Eq. (2) is 1. The last layer gathers the hidden layer's outputs and readjusts 10 the output to provide the correct function value. A typical input would activate several neurons (the weight is not exactly the same as the input since the input is not in the training set), i.e. 0<a<1 for several neurons output and the final output to the network is a combination of the different neuron outputs. This feature allows to the network to interpolate the function values and, therefore learn the shape of the function. Assum- 15 ing that the training set is well spread along the input range, the RBF network learns the shape of the function with the training set. An independent set, the test set, is used to evaluate the function approximation. It works like a spline interpolation with the advantage that N-dimensional functions can be approximated easily but the disadvantage that the function is unknown.

Development of a neural network-based model for the AOD estimation
We need to relate the wavelengths of the WSI filters with the wavelengths used to measure the direct irradiance with the CIMEL CE318. We used the nearest wavelength and, therefore the 450 nm filter is associated with the 440 nm in CIMEL, the 650 nm filter is associated with the 675 nm in CIMEL and the 800 nm filter is associated with the 25 870 nm in CIMEL. Thus, we developed a RBF network model for each WSI wavelength obtaining an estimation of the AOD at 440, 675 and 870 nm. 19998 Inputs to the RBF network are the SZA • corresponding to the measurement and the radiance at several scattering angles over the principal plane.

Neural network-based model
The whole data set consists of the radiances over the principal plane of the 1047 images, the SZA • at the time of the measurement and the synchronous CIMEL measure-5 ment. This set is divided randomly in two subsets, a training set using 2/3 of the whole data set and a test set using 1/3 of the whole data set. All values are normalized to the range [0,1], i.e. the values are rescaled to that range where the minimum and maximum values correspond to 0 and 1 respectively. Inputs (radiances over the principal plane and SZA • ) and outputs (AODs) are normalized. During the training every measurement is used to adjust the internal values of the neural network (the weights) so the output is the desired variable amount, i.e. the AOD corresponding to the wavelength we are trying to estimate. Once the network is trained we use the test set to evaluate the performance of the network. The coefficient of determination (R 2 ) calculated between the AOD measured with the CIMEL and the values estimated with the network using the 15 test set is our estimator of the performance. The performance depends on the selection of the training and test sets. These sets are created randomly out of the whole data set. Hence we repeat the process nine times and select the best partition, i.e. the one that yields the best performance. Interactive Discussion principal plane, but not all the scattering angles are necessary. The first iteration of the algorithm creates 100 networks using the radiance at one of the scattering angles. All of them are evaluated and the best one is added to the solution. The second iteration of the algorithm creates 99 networks using the radiance at the best scattering angle of the previous iteration and the radiance at a different scattering angle. Once again, 5 the best one is added to the solution. This process is repeated until the performance decreases, i.e. no more scattering angles are needed. During the process we applied a mask to eliminate the measurements with non valid values (obstruction due to the horizon or shadow system). As a result, the final data set may vary depending on the scattering angles used.

Results
Every wavelength has an independent model. The inputs and, therefore the training and test sets, are different for each AOD estimation.
The greedy algorithm selected only one scattering angle for each AOD model. For the blue WSI channel (450 nm associated to the 440 nm AOD), it selected the 37 • 15 scattering angle. Therefore the model has two inputs: the radiance of the sky at that scattering angle over the principal plane and the SZA • . 117 measurements had to be eliminated from the original 1047 measurements due to shadow system obstruction, so the model was created from the remainder of this set with 930 measurements. The greedy algorithm for the model with the red WSI channel (650 nm associated 20 to the 675 nm AOD) selected the 71 • scattering angle. The number of valid measurements after applying the mask is 968. The greedy algorithm for the model with the NIR WSI channel (800 nm associated to the 870 nm AOD) selected the 83 • scattering angle. The number of valid measurements is 973 in this case. provides information about the over -or underestimation associated with the model. The coefficient of determination provides an evaluation of the experimental variance explained by the model. The root mean square deviation (RMSD) and the mean bias deviation (MBD) have been also evaluated as quality estimators. These quality estimators allow us to evaluate the differences between the experimental data and the 5 model and the presence of a systematic over -or underestimation. Table 1 shows the statistics for the three models. Figure 7 show histograms of the differences between measured and estimated AOD at 440, 675 and 870 nm.
As we can see in Fig. 6a and Table 1, 96% of the data variance is explained by the model that estimates the AOD at 440 nm. MBD and the slope of the linear fit reveal 10 a slight systematic underestimation. Figure 7 shows a histogram of the differences between the calculated and estimated values. It reveals that 81% of the estimated AOD values at 440 nm had a deviation less than 0.01 with respect to the CIMEL result which is the AERONET AOD estimated uncertainty (Holben et al., 1998). Figure 6b and Table 1 reveal that 94% of the data variance is explained by the model 15 that estimates the AOD at 675 nm. MBD and the slope of the linear fit also indicate a slight systematic underestimation. Figure 7 shows that almost 80% of the estimated AOD at 675 nm has a deviation less than 0.01. Finally, Fig. 6c and Table 1 reveal that 92% of the data variance is explained by the model that estimates the AOD at 870 nm. While MBD suggests a slight overestimation, 20 the slope of the linear fit indicates an underestimation. Figure 7 shows that 90% of the estimated AOD at 870 nm has a deviation less than 0.01. Figure 6a, b and c reveal that the data set is not homogeneously distributed along the whole range of values. There are a lot of points with low AOD and very few with higher AOD. This can explain the slight systematic underestimation of the model. The 25 linear fit is forced to zero and there are a lot of points close to zero, but the very few values far from the zero introduce a variance that, in this case makes the slope be slightly below 1.
The coefficient of determination decreases when we estimate AOD at longer wave-20001 works as seen in Sect. 3.3, and β parameter is also estimated. Secondly, a new neural network has been trained using the calculated α by AERONET. The AERONET algorithm calculates α using a different interval of wavelengths. The interval 440-870 nm includes the values we estimate, and therefore this is the interval used to compare the results in both approaches. 15 For the first method, Fig. 8a shows estimated versus calculated values of α with the CIMEL in the interval 440-870 nm using the standard AERONET procedure. Figure 8c shows the histogram of the differences between calculated and estimated α. Figure 8a reveals that 63% of the data are explained by the model that estimates α. Figure 8c shows that 48% of the estimated α has a deviation less than 0.1, which is 20 the estimated uncertainty in the AERONET procedure for α calculation (Holben et al., 1998).
The α estimation is affected by the error introduced in the AOD estimation. The AOD at 870 nm introduces an error in the calculation of α by linear fitting of ln(AOD) vs. wavelength. For this reason, we tried a new neural network-based model using RBFs 25 to estimate the value of α using the radiance over the principal plane of the sky images for the three wavelengths. The inputs for this model are the same for the estimation of the AOD together, i.e. we combined the radiance at different scattering angles over 20002 data are explained by the model that estimates α. Figure 8c shows that 84% of the estimated α has a deviation less than 0.1, which is the estimated uncertainty in the AERONET procedure for α calculation. This represents a clear improvement in the estimation of α from WSI images. Even though the uncertainty in the estimation of α is large with the standard method, 10 the estimation is still useful for the interpolation of the AOD at different wavelengths.
We have tested it calculating the AOD at 500 nm with the α and β estimated with the first method usingÅngström Law (Ångström, 1964) and compared it with the AOD at 500 nm obtained from CIMEL CE318 measurements. Figure 9 shows estimated values of AOD at 500 nm versus calculated values of AOD with the CIMEL. 96% of the data 15 variance is explained by the model that estimates the AOD at 500 nm. As we can see, the use of this estimation yields a very good estimation of AOD in different wavelengths (one of the main uses of α). However, if we need a more precise estimation of α, for example as input in an inversion model, we also have the most accurate estimation using the neural network-based model.

Conclusions
Three neural network-based models using RBFs has been created to estimate the value of the AOD at three different wavelengths using the radiance over the principal plane of the sky images from the calibrated sky imager WSI. The correlation constants are close to unity and the number of cases within measurement error is very large.
performed in two ways. First, it has been calculated using the standard procedure in AERONET using the AODs estimated with the neural network (the β parameter is also calculated) and secondly, it has been estimated with a new neural network-based model using RBFs. Inputs to this RBF are the radiances in the same scattering angles used for the AOD models.

5
The three AOD models provide an estimation that, according to validation, is inside the nominal error of the AERONET (±0.01) in approximately 80% for the blue and red channels and 90% for the NIR channel. In all cases the models explain up to 92% of the variance of the experimental data. The coefficient of determination decreases when we estimate AOD at longer wavelengths. This can be caused by the difference between the 10 central wavelength of the filter in the sky imager and the CIMEL. This difference is larger with longer wavelengths and, therefore, the overlap of the wavelength decreases and so the performance. Nevertheless, all estimations have a coefficient of determination over 0.92.
Concerning the scattering angles selected with the greedy algorithm, these reveal 15 that to estimate the AOD at 450 nm (blue) it is necessary to use a point close to the sun. The estimation of the AOD at 870 nm (NIR) requires a point farther from the sun. At 675 nm the behavior can be catalogued as in the middle of the behavior presented in the other two wavelengths. The α estimation using the Aeronet method is affected by the cumulative error of all 20 the AOD estimations. However, almost 50% of the data are inside the nominal error of the AERONET program for α calculation using the standard procedure in AERONET. The neural network-based model for α estimation increases the explanation of the data to 63% and the data inside the nominal error is increased to 77%. The neural network process is complex but increases substantially the estimation. The neural network 25 model forces the result to be the same as that estimation of α for the interval 440-870 nm.
The results are promising in the sense that it seems to be feasible that a sky imager can estimate the AOD and the algorithm could be applied in the field. Furthermore, the 20004 Government and his research stay at University of California at San Diego has been also funded by the Andalusian Regional Government. We are also especially thankful to the principal investigator of the AERONET site at SGP, Rick Wagener and the collaboration of the US Department of Energy as part of the Atmospheric Radiation Measurement Program Climate Research Facility.
properties assessed by an inversion method using the solar principal plane for non-spherical particles, Journal of Quantitative Spectroscopy and Radiative Transfer, 109(8), 1504-1516, 2008. 20007 Table 1. Statistical results of the validation for the radial basis networks for estimation of AOD at 440 nm, AOD at 670 nm and AOD at 870 nm. The column B represents the slope of the linear fitting through zero of the data, R 2 is the coefficient of determination, MBD is the mean bias deviation and RMSD is the root mean squared deviation.