CLOUD PROBABILITY MASK FOR PROBA/CHRIS HYPERSPECTRAL IMAGES
Luis G´ omez-Chova1 , Gustavo Camps-Valls1 , Julia Amor´ os-L´ opez1 , Javier Calpe1 , Luis Guanter2 , Luis Alonso2 , Juan C. Fortea2 , and Jos´ e Moreno2 2
1 GPDS, Electronic Engineering Department, University of Valencia. LEO, Earth Sciences and Thermodynamics Department, University of Valencia. Doctor Moliner 50, 46100, Burjassot (Valencia), Spain. Corresponding author e-mail:
[email protected]
Abstract This paper presents a methodology that faces the problem of detecting clouds in multispectral images acquired by sensors working in the VNIR range. The method is capable of detecting clouds accurately providing probability and cloud abundance rather than flags. For this purpose, the method exploits the special characteristics of the PROBA/CHRIS instrument, such as the high number of spectral bands, the spatial resolution, and the specific band locations. Performance of the proposed algorithm is tested on six CHRIS images, which represent scenarios with different characteristics. Results demonstrate that the proposed algorithm yields good results even classifying difficult cloud pixels, especially thin cirrus clouds and clouds over ice/snow, and provides a better description of detected clouds. Key words: cloud detection, cloud masking, hyperspectral images, PROBA, CHRIS, unsupervised classification, spectral unmixing.
1.
INTRODUCTION
In remote sensing images acquired by instruments working in the visible and near-infrared range (VNIR) of the electromagnetic spectrum, undetected clouds are the most significant source of error for true ground reflectance estimation [1]. Therefore, the operational use of these images can be hampered by the presence of clouds. Cloud screening approaches heavily depend on the characteristics of each sensor. Obviously, its spectral and spatial resolution, or the spectral range are critical. In particular, the Compact High Resolution Imaging Spectrometer (CHRIS) [2] provides us with a number of channels in the visible (VIS) and near infrared (NIR) range with enough spectral resolution to increase the cloud detection accuracy [3]: it covers a spectral range from 400 nm to 1050 nm with a maximum spatial resolution of 17 or 34 m at nadir depending on the acquisition mode (more information on http://earth.esa.int/proba/). This high resolution makes CHRIS a good choice in order to propose and validate cloud detection methodologies. CHRIS is mounted on board the European Space Agency (ESA) small satellite platform called PROBA (Project for On Board Autonomy). Thanks to the PROBA platform pointing capabilities, the acquisition plan tries to avoid acquisitions with cloud coverage, but occasionally images are partially affected by clouds. In these cases, users that requested the acquisition have a special interest on both cloud location and contribution to the registered spectra. In this context, the main objective of this work is to develop and validate a method for cloud detection using self-contained information provided by PROBA/CHRIS data. The method should be capable of: detecting clouds accurately; providing probability and cloud abundance instead of flags; better describing detected clouds (cloud abundance, type, height, subpixel coverage); and providing information about shadowed areas in order to take into account the spectral changes. As a consequence, by masking only the image areas affected by cloud covers, the whole image is not necessarily discarded, thus making multitemporal and multiangular studies possible.
Figure 1. RGB composite of the CHRIS images (Mode 1, FZA 0o ) over the test sites of BR-2003-07-14, BR2004-05-27, BR-2005-07-17, OH-2005-01-08, LG-2004-02-29, and RC-2004-04-23 (images displayed in rows from left to right). 2.
DATA MATERIAL
For this study, a dataset consisting of six CHRIS Mode 1 acquisitions (five multi-angular images per acquisition) over four of the core sites of the PROBA mission was considered. In particular, images taken over Barrax (BR, Spain), Oahu (OH, Hawaii, USA), Lake Argyle (LG, Australia), and Reynold’s Creek (RC, USA) have been included in the study (Fig.1) in order to take into account their different characteristics: geographic location (latitude/longitude); date and season; type of cloud (cumulus, cirrus, stratocumulus); and surface types (soil, vegetation, sand, ice, snow, lakes, sea, etc). These images are useful to validate the performance of the algorithm since they include different landscapes and two critical cases given the especial characteristics of the induced problems, ice and snow.
2.1.
IMAGE PRE-PROCESSING
All the images used in this work have been reprocessed by SIRA to the version 4.1 of the CHRIS HDF files [4], which presents an improved radiometric calibration [5]. In addition, the information contained in the metadata attributes of the CHRIS HDF file (image date, azimuth and zenith angles, etc.) and the quality mask added to the HDF data (pixel saturation and occurrence of errors) are used during the image pre-processing. One of the errors affecting CHRIS images is the fact that transmission of channel 2 randomly produces wrong odd pixels in image rows. We correct this drop-out noise by substituting the wrong pixel value by a weighted spatio-spectral average of its nearest neighbors. In order to avoid the poor performance of spatial filters in border or inhomogeneous areas, we weight the contribution of each neighboring pixel by its spectral distance to the corrected pixel. It is worthnoting that the values of bands with errors (indicated by the CHRIS HDF image quality mask) are not included in the computation of the spectral distance (Euclidean distance). Another well-known problem is a spatial coherent noise known as vertical stripping (usually found in images acquired by push-broom sensors). This multiplicative noise in image columns comes from irregularities of the entrance slit and CCD elements in the across-track direction. In [6], the authors proposed a method to remove and model the vertical stripping of CHRIS instrument. The main advantage of this method is that it can distinguish between high spatial frequencies coming from the surface contribution and the vertical stripping. Finally, CHRIS products are provided in top of the atmosphere (TOA) radiance (radiometrically calibrated data). Therefore, TOA reflectance is estimated to reduce the data dependence on particular illumination
CHRIS Mode1 spectral channels 1
CHRIS Mode1spectral bands Atmospheric transmission Vegetation spectral signature Bare soil spectral signature
0.5
0 400
600 800 wavelength (nm)
1000
Figure 2. CHRIS Mode 1 band locations (red) superimposed to the atmospheric transmittance (blue) and a reflectance spectrum of healthy vegetation (green) and bare soil (black). conditions (day of the year and angular configuration) and illumination effects due to rough terrain. TOA apparent reflectance is estimated according to: ρ(x, y, λ) =
πL(x, y, λ) , cos(θI )I(λ)
(1)
where L(x, y, λ) is the provided at sensor upward TOA radiance, I(λ) is the extraterrestrial instantaneous solar irradiance, and θI is the angle between the illumination direction and the vector perpendicular to the surface. In this work, θI is approximated by the Solar Zenith Angle provided in the CHRIS HDF file attributes since one can assume flat landscape and constant illumination angle for the area observed in a CHRIS image.
3.
CLOUD SCREENING ALGORITHM
The standard approach to detect clouds in a scene is the use of a set of static or adaptive thresholds [7] applied to some features (brightness, temperature, etc.) in every image pixel. However, these methods can fail for several reasons such as subpixel clouds, land covers with high reflectance, illumination and observation geometry, sensor calibration, variation of the spectral response of clouds with cloud type and height, etc. [1]. In this work, we present a cloud screening procedure which is presently under refinement and extensive validation [8]. Essentially, the cloud detection algorithm is constituted by the following steps.
3.1.
Feature Extraction
The measured spectral signature depends on the illumination, the atmosphere, and the surface. Fig.2 shows the location of CHRIS Mode 1 bands compared with the atmospheric transmittance, and spectral signature of healthy vegetation and bare soil. Physically-inspired features that increase separability of clouds and surface covers can be extracted independently from the bands that are free from atmospheric effects and from the bands affected by the atmosphere. For illustration purposes, the figures of this section will show maps of the extracted features from the CHRIS-BR-040527-416E-41 image. Regarding the reflectance of the surface, one of the main characteristics of clouds is that they present bright and white spectra (Fig.3). • A bright spectrum means that the intensity of the spectral curve (related to the albedo) should present R relatively high values. Therefore, cloud brightness is calculated as the integral of spectrum, ρ(λ)dλ. • A white spectrum means that the first derivative of the spectral curve should present low values. For the present method, the average of the spectral derivative, hdρ(λ)/dλi, has been chosen as one of the features. Regarding the atmospheric absorptions present in the spectrum of a pixel, another meaningful feature is the fact that clouds are at a higher altitude than the surface. It is worth noting that atmospheric absorption depends on the atmospheric components and the optical path. Since light reflected on high clouds crosses a shorter section of the atmosphere, the consequence would be an abnormally short optical path, thus weaker atmospheric absorption features. Atmospheric oxygen or water vapour absorptions (at 760 nm and 940 nm respectively) can be used to estimate this optical path.
Spectral intensity
Mean spectral derivative
−3
x 10
0.25
1
0.2 0.5
0.15 0
0.1 −0.5
0.05 −1
Figure 3. Cloud brightness (left) and the mean spectral derivative (right) of the TOA reflectance of the CHRISBR-040527-416E-41 image. Normalized Atmospheric Transmittance in the Oxygen absorption band
Oxygen absorption estimation 0.55
1
0.5
0.45
0.5
755
760
765
770
775
Spectral efficiency of CHRIS channels 1.2
TOA reflectance
0.4
0
0.35
0.3
0.25
1
0.2
0.8
0.15
0.6
0.1
0.4 755
760
765 wavelength (nm)
770
775
0.05 400
500
600
700 wavelength (nm)
800
900
1000
Figure 4. Intermediate products used to estimate the optical path, d, from the oxygen absorption band. Left: effective atmospheric vertical transmittance, exp(−τ (λ)), estimated for the CHRIS channels. Right: the interpolated, ρ0 (λ), and measured, ρ(λ), reflectance inside the oxygen band, and the estimated reflectance at the maximum absorption (λ=760.625 nm). The light transmitted through a non-dispersive medium can be expressed by (Bouguer-Lambert-Beer law): L(λ) = L0 (λ) exp(−τ /µ) = L0 (λ) exp(−τ (λ)d/µ) ,
(2)
where L0 (λ) is the initial light, the term exp(−τ /µ) is the transmittance factor, 1/µ is the optical mass, τ is the optical depth, τ (λ) is the atmospheric optical depth for a vertical path, and d is a factor accounting for the path crossed by the radiation (d can be interpreted as a product of the component concentration and the distance crossed by the radiation, which will be approximately 1 when the light crosses one atmosphere with a vertical path). Assuming that there are no horizontal variations in the atmospheric component concentrations, one can consider that the optical path only depends on the height. Therefore, an estimation of the atmospheric absorption provides a measure of the optical path: d(x, y) = µ
− ln(L(x, y, λi )/L0 (x, y, λi )) , τ (λi )
(3)
where the effective atmospheric vertical transmittance, exp(−τ (λ)), is estimated from a high resolution curve for the channels λi of the instrument affected by the absorption; and the spectral signature without absorption, L0 (x, y, λi ), is estimated by interpolating the nearby channels that are unaffected by this absorption. The approach followed in this paper can be devised from Fig.4 for the so-called Oxygen-A band. The foundations for the procedure are solid but the performance of the instruments severely affects the robustness of the method. In the case of the oxygen absorption, at least one CHRIS Mode1 spectral band is affected by the absorption. However, its wavelength is not necessarily centered in the maximum absorption and its bandwidth is too high making difficult an accurate estimation of the optical path. Moreover, residual vertical striping noise is critical when ratios with low values are computed in (3). Therefore, when using the oxygen absorption band (760 nm), only high clouds are well distinguished from the surface and low clouds present almost the same values as the surface. In the case of the water vapour absorption (maximum absorption at 940 nm close to the end of CHRIS spectral range), the water vapour distribution is extremely variable. However, one can assume small differences since CHRIS spatial coverage is small (15 km swath), and almost all the atmospheric water vapour is distributed in the first 2-3 km of the atmosphere below most of the cloud types. Therefore, the optical path estimated from this absorption provides better cloud discrimination since this band is broader than the oxygen band and more CHRIS bands are completely affected (Fig.5).
Oxygen absorption
Water Vapour absorption 2
0.3
1.8 1.6
0.25
1.4
0.2 1.2 1
0.15
0.8 0.6
0.1
0.4
0.05 0.2 0
Figure 5. Estimation of the optical path from the oxygen (left) and the water vapour (right) absorption band for the CHRIS-BR-040527-416E-41 image. 3.2.
Region of interest
The set of representative features for the cloud detection problem is extracted in order to improve the segmentation of the image. In addition, results of the clustering algorithm will also improve when applied over the regions of the image where clouds are statistically representative. In order to find regions that potentially could contain clouds, hard non-restrictive thresholds are used to provide a first map of cloud-like pixels. Then, a region growing algorithm is carried out, along with a morphological process that dilates cloudy areas. The result is not intended to be a classification, but to obtain a mask or region of interest (ROI), in which presence of clouds is significant for the later clustering.
3.3.
Image Clustering and Labeling
The clustering algorithm is applied only to the selected ROI. The inputs for each pixel, xk , are the extracted features: the brightness in the VIS and the NIR region, the mean spectral derivative of the TOA reflectance, and the estimated oxygen and water vapour absorptions. We use the Expectation-Maximization (EM) algorithm [9] to estimate the parameters of a Gaussian mixture model since it considers the full relationship among variables and provides probability maps for each cluster. If the ROI has been correctly found, even using a low number of clusters, N , some of them should correspond to different cloud types. Once clusters are determined, the spectral signature of each cluster, si (λ), is estimated as the average of the spectra of the cluster pixels. At this point of the process, the obtained clusters can be labeled into geo-physical classes taking into account three complementary sources of information (Fig.6): the thematic map with the spatial distribution of the clusters in the scene, the spectral signatures of the clusters, si , and the location in the image of the pixels with the spectral signature closer to si . This information can be either analyzed directly by the user or compared to a spectral library with representative spectra of all the classes of interest. Once all clusters have been related to a class with a geo-physical meaning (Fig.7, left), or at least they have been classified as cloud or non-cloud, it is straightforward to merge all the clusters belonging to a cloud type. Since PN the EM algorithm provides posterior probabilities (Pik ∈ [0, 1] and i=1 Pik = 1), a probabilistic cloud index, based on the clustering of the extracted features, can be computed as the sum of the posteriors of the cloudP clusters: Cloud Probability k = i Pik for all cluster i classified as cloud (Fig.7, center ). However, if clusters are well separated in the feature space, the posteriors decrease drastically from one to zero in the boundaries between clusters (Fig.7, center ). Therefore, this Cloud Probability index indicates when one pixel belongs more to a cloud-cluster than to one of the other clusters found in the image, but it does not give information about the cloud content at subpixel level (crucial when dealing with thin clouds or partially covered pixels).
3.4.
Spectral Unmixing
In order to obtain a cloud abundance map for every pixel in the image, rather than flags or a binary classification, a spectral unmixing algorithm is applied to the multispectral image using the full spectral information. In our case, we are going to consider the spectral signatures of the clusters, si , as the representative pixels of the covers present in the scene (endmembers) , and they are used to build matrix M. The Fully Constrained Linear Spectral Unmixing (FCLSU) [10] algorithm is used to perform the spectral unmixing where each pixel-k is
Cluster centers
Clustering of selected areas 1.4
1 2 3 4 5 6 7 8
8 1.2
6 5 4 3
TOA Apparent Reflectance
7 1
3
1
4 7
0.8
8 2
0.6
0.4
2
6 5
0.2
1 0
0 400
500
600
700 800 wavelength (nm)
900
1000
1100
Figure 6. Thematic map with the distribution of the clusters in the scene (left), the spectral signatures of the clusters (center), and the location in the image of the pixels with the most similar spectra (right). Cloud Probability
Classification of selected areas
Unmixing Cloud Abundances 1
1
0.9
0.9
0.8
0.8
0.7
0.7
5:clouds
0.6
0.6
4:shadows
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
8:shadows
7:soil
6:bright clouds
3:soil
2:vegetation
1:cirrus
0:background
Figure 7. Thematic map with the distribution in the scene of the classes of the clusters (left), Cloud Probability index computed from the posteriors of the cloud-clusters (center); and Cloud Abundance computed from the unmixing coefficients of the cloud clusters (right). modeled as a mixture of the cluster centers: ρk = M · ak + ε. The vector ak contains, for the sample pixel k, the abundances of the spectral signatures of the clusters, which are related to a class with a geo-physical meaning. As it happens with the probabilities of the clusters, the abundance values are constrained to aik ∈ [0, 1] and PN of cloud is computed as the sum of the abundances of the cloudi=1 aik = 1. Therefore, the abundance P clusters: Cloud Abundance k = i aik for all cluster i classified as cloud (Fig.7, right). An improved cloud abundance map can be obtained when combining the Cloud Abundance and the Cloud Probability by means of a pixel-by-pixel multiplication. That is, combining two complementary sources of information processed by independent methods: the cloud fraction or mixing (obtained from the spectra) and the cloud probability that is close to one in the cloud-like pixels and close to zero in remaining areas (obtained from the extracted features).
4.
RESULTS
We have tested the proposed scheme for the aforementioned images. All the method has been implemented in Matlab programming language and it is capable to handle directly the CHRIS HDF original files and process them automatically producing the results in an ENVI standard format. In the previous section, intermediate images of the method have been shown (Fig.3 to Fig.7). This section analyzes the performance of the method on the six images by means of the thematic map of the clusters and the cloud abundance product (Fig.8). The following conclusions about the method results are depicted in Fig.8. The three images over the Barrax (BR, Spain) site are a good example of an easy cloud detection problem, when clear clouds are well contrasted with soil and vegetation. The ROI selection can be easily appreciated in the classification images, being more important in the BR-2003-07-14 image where small clouds could be mixed in a cluster with other classes if the whole image is considered. At the OH-2005-01-08 image (Oahu, Hawaii), waves and beach sand pixels belong to a cluster labeled as cloud due to their high reflectance and whiteness, but they present low probabilities and abundances. In addition, this image presents cirrus clouds over land and over sea that are well detected since specific clusters describe them and have been correctly classified. One of the weak points of the algorithm is the use of thresholds to select the ROI, because some thin or small clouds can be excluded from the ROI. A solution is to relax thresholds at the risk of considering the whole image as ROI as in the image of LG-2004-02-29 (Lake Argyle, Australia). But, even in this case, results are good if clouds cover a sufficient percentage of the image
Figure 8. Classification (first and third rows) and cloud product (second and fourth rows) for BR-2003-07-14, BR-2004-05-27, BR-2005-07-17, OH-2005-01-08, LG-2004-02-29, and RC-2004-04-23 (from left to right). or the number of clusters is high enough. The image of RC-2004-04-23 (Reynold’s Creek, United States) is an example of one of the critical issues in cloud detection, the presence of ice/snow in the surface. These pixels and clouds have a similar reflectance behavior. However, the atmospheric absorption suffered by cloud pixels is lower than for the surface pixels due to the their height, thus different clusters are found for these two classes in the image. Thanks to the extracted atmospheric features, ice/snow pixels present low Cloud Probability values although the Cloud Abundance provided by the spectral unmixing could be relatively high due to the spectral similarities. In consequence, both information types are combined improving the final classification accuracy.
5.
CONCLUSIONS
In this paper, we have presented a new technique that faces the problem of accurately identifying the location and abundance of clouds in hyperspectral images. At present, we are refining the proposed scheme and testing its performance on PROBA/CHRIS data. In any case, the procedure should serve as a test to develop a cloud screening algorithm for other spectral sensors in the range from 400 nm to 1000 nm. The use of this data allow us to assess algorithm performance in favorable spatial resolution (34 m) and number of bands (62 channels). Further refinements should be introduced in order to enhance the robustness of the procedure. For instance, dynamic thresholds should be preferred to find the regions to analyze; final maps can be further analyzed through texture algorithms; and edge detection of the clouds and shadow regions, and then apply morphology operators that extract information on the structure of the image. It has been observed that co-registration and geo-coding of images are critical in order to introduce a series of improvements: (1) use DEM to correct the image and estimate heights, (2) combine the information provided by the five images thus including multi-angularity, (3) take into account sun position in order to relate cloud and shadow positions, and (4) use a coastline map in order to considerate separately clouds over land or water.
ACKNOWLEDGMENTS This paper has been partially supported by the Spanish Ministry for Education and Science under project ESP2005-07724-C05-03 and by the Generalitat Valenciana under project HYPERCLASS/GV05/011. Authors wish to thank ESA and SIRA for the availability of the image database and the assistance provided.
REFERENCES 1. Simpson, J. Improved cloud detection and cross-calibration of ATSR, MODIS and MERIS data. In ESA-SP479, editor, ATSR International Workshop on the Applications of the ERS along track scanning radiometer, ESRIN, Frascati, Italy, June 1999. 2. Barnsley, M., Settle, J., Cutter, M., Lobb, D., and Teston, F. The PROBA/CHRIS mission: a low-cost smallsat for hyperspectral, multi-angle, observations of the Earth surface and atmosphere. IEEE Trans. Geosci. Remote Sensing, 42(7):1512–1520, July 2004. 3. Diner, D., Clothiaux, E., and Girolamo, L. D. MISR Multi-angle Imaging Spectro-Radiometer Algorithm Theoretical Basis. Level 1 Cloud Detection. Jet Propulsion Laboratory, JPL D-13397, 1999. 4. Cutter, M. and Johns, L. CHRIS data products – latest issue. In ESA-SP-593, editor, 3rd CHRIS/Proba Workshop, ESRIN, Frascati, Italy, June 2005. 5. Cutter, M. CHRIS calibration issues. In ESA-SP-578, editor, 2nd CHRIS/Proba Workshop, ESRIN, Frascati, Italy, April 2004. 6. G´omez-Chova, L., Alonso, L., Guanter, L., Camps-Valls, G., Calpe, J., and Moreno, J. Modelling spatial and spectral systematic noise patterns on CHRIS/PROBA hyperspectral data. In Bruzzone, L., editor, Image and Signal Processing for Remote Sensing XII, volume 6365, page 63650Z. SPIE, 2006. 7. Vittorio, A. D. and Emery, W. An automated, dynamic threshold cloud-masking algorithm for daytime AVHRR images over land. IEEE Trans. Geosci. Remote Sensing, 40(8), 2002. 8. G´omez-Chova, L., Amor´ os, J., Camps-Valls, G., Mart´ın, J., Calpe, J., Alonso, L., Guanter, L., Fortea, J., and Moreno, J. Cloud detection for CHRIS/Proba hyperspectral images. In Schafer, K. and Comeron, A., editors, Remote Sensing of Clouds and the Atmosphere X, volume 5979, page 59791Q. SPIE, 2005. 9. Dempster, A. P., Laird, N. M., and Rubin, D. B. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B, 39:1–38, 1977. 10. Chang, C. Hyperspectral Imaging: Techniques for Spectral Detection and Classification. Kluwer Academic/Plenum Publishers, 2003.
