CIRRUS REMOVAL IN MULTISPECTRAL DATA WITHOUT 1 . 38 μ M SPECTRAL DATA

Cirrus is one of the most common artifacts in the remotely sensed optical data. Contrary to the low altitude (1-3 km) cloud the cirrus cloud (8-20 km) is semitransparent and the extinction (cirrus influence) of the upward reflected solar radiance can be compensated. The widely employed and almost ’de-facto’ method for cirrus compensation is based on the 1.38 μm spectral channel measuring the upwelling radiance reflected by the cirrus cloud. The knowledge on the cirrus spatial distribution allows to estimate the per spectral channel cirrus attenuation and to compensate the spectral channels. A wide range of existing and expected sensors have no 1.38 μm spectral channel. These sensors data can be corrected by the recently developed haze/cirrus removal method. The additive model of the estimated cirrus thickness map (CTM) is applicable for cirrus-conditioned extinction compensation. Numeric and statistic evaluation of the CTM-based cirrus removal on more than 80 Landsat-8 OLI and 30 Sentinel-2 scenes demonstrates a close agreement with the 1.38μm channel based cirrus removal.


INTRODUCTION
The mean estimate of the global cirrus distribution is 16.7% (Sassen et al., 2008) and the maximum value obtained is over 50% (Heidinger andPavolonis, 2005, Chepfer et al., 2000).Such a high expectation rate of cirrus in the atmosphere defines a relatively high probability of cirrus in remotely sensed optical data.The radiance extinction conditioned by cirrus has an important property: the surface reflected solar radiance is partly attenuated leaving an opportunity to be restored.The 1.38µm spectral channel in a multispectral instrument design gives an opportunity to estimate the thickness distribution of cirrus cloud and to restore the upwelling reflected radiance for semitransparent cirrus.The state-of-the art methods for cirrus correction employ the 1.38µm spectral channel for the restoration (B.C. Gao et al., 2002, Gao et al., 2004, Richter et al., 2011).
The absence of the 1.38µm cirrus channel in a huge set of existing and planned multi/hyperspectral sensors requires a correction method independent of the cirrus channel.In this paper we demonstrate cirrus removal by the method for haze/cirrus compensation (Makarau et al., 2014).Instead of the haze thickness map (HTM) the Cirrus Thickness Map (CTM) can be calculated for semi-transparent cirrus cloud and spatially varying cirrus thickness.The spatial and spectral cirrus distribution is estimated and the multispectral data are restored with the use of the additive model.Cirrus map and cirrus corrected (decirrus) data are employed as the reference.The accuracy of cirrus correction is numerically and statistically evaluated.

METHOD
We use an additive model to calculate the cirrus influence: where L sensor is the acquired radiance, L0 is the sum of path radiance and surface reflected radiance, and CR is the cirrus contribution.In terms of a sensor acquired digital number (DN) (a linear model of the DN to radiance conversion) Eq.1 can be written as where (x, y) are the pixel coordinates in an image, i is the spectral channel number, DN sensor i (x, y) is the acquired DN, and DNi(x, y) is the DN without the attenuation by cirrus.The term CR ′ i (x, y) is dependent on the cirrus thickness and defined as cirrus thickness map CT Mi(x, y).In this paper the cirrus thickness is defined in the range of the spectral channel radiance or the DN.The subtraction of the CT M allows to restore the data.

Spatial Cirrus Thickness
The CT M thickness is estimated by searching for (non-overlapping window w × w) dark pixels in a channel with a minimal ground reflectance and a maximal attenuation by cirrus.A blue band in the region 0.37-0.49µm is most suitable.The CT M estimate given a single channel leads to an overcorrection of this channel.This is overcome by a new synthetic spectral channel creation (linear extrapolation is used): Bex(x, y) = EXT RAP OL(Bm(x, y), Bn(x, y)), (3) where Bex(x, y) is the extrapolated spectral channel, Bm(x, y) and Bn(x, y) are two blue/green spectral channels from the data.The CT M (x, y) estimate given a lower wavelength Bex is more precise since surface reflectance is less and cirrus thickness is higher.Large surface areas with a high albedo make the selection of dark pixels difficult using a small window.This area is labeled by a segmentation method and excluded from the search.

Spectral Cirrus Thickness
Cirrus thickness in terms of sensor measured radiance usually decreases from the shortest to the longest wavelength.Subtraction of a constant CT M (x, y) from all the channels leads to a loss of spectral consistency.A channel-specific CT Mi(x, y) is calculated and the slope (ki) of the linear regression with the dependent CT Mi(x, y) versus the independent CT M (x, y) (calculated using the Bex) is stored in ki.The ki coefficient is used to scale the CTM.The regression is calculated using the pixel values in the cirrus regions: where ki is the slope of the fitted line for band i, K is the array of the coefficients (ki ∈ K).

Cirrus correction and Clear Scene Aerosol Compensation
A subtraction of the CT M (x, y) * ki from the band DN sensor i (x, y) recovers the decirrus band: The subtraction removes the cirrus along with possibly a fraction of the clear scene aerosol thickness and the data should be compensated for the clear scene aerosol fraction.CT M (x, y) pixels outside cirrus contain values of the clear scene aerosol thickness (the aerosol thickness is also calculated in cirrus-free areas).A compensation for the aerosol thickness is performed: where the cirrus and no cirrus are cirrus and non cirrus pixels, respectively.

NUMERIC EVALUATION
The data restored by the 1.38µm band decirrus method is employed as the reference.A small error in this restoration always exists originating from the methodology and the extincted radiance which can never be restored with 100% accuracy.The same error assumptions are valid for the proposed CTM method.The experiments were carried on a database of Landsat-8 and Sentinel-2 scenes, in total more than 120 scenes.We demonstrate restoration results for one Landsat-8 and one Sentinel-2 scenes.
The numeric evaluation is expected to show a small difference between the 1.38µm band decirrus method and the CTM to prove the hypothesis that the CTM restores the data with a comparable accuracy.The evaluation employs the absolute difference histogram (the histogram of the 1.38µm method and the CTM decirrus results difference), the relative difference histogram (the histogram of the two decirrus results difference divided by the reference result), the restored data histograms, and spectrum profile taken in a pixel in a selected cirrus restored area.Cirrus maps and the RGB compositions (subscenes) are presented for a visual interpretation.1).The relative difference peak value is 6% demonstrating a small difference for this measure (Figure 2 c) and Table 1).The spectra collected in the cirrus area (Figure 3) demonstrate spectral consistency for the 1.38µm and the CTM decirrus methods.The CTM method compensates slightly less cirrus in the 1.6 and 2.2µm band.This is an expected behaviour since the extinction in these bands is primarily conditioned by absorption and less by scattering.Generally, the 1.38µm based decirrus method performs with slightly less accuracy in the atmospheric absorption regions (e.g. in the 0.94 and 1.38µm) compared to the visible and near infrared (Richter et al., 2011) and these bands are omitted.Mean values of the absolute and relative difference as well as the root mean squared error (RMSE) calculated for each band are presented in Table 1 demonstrating that only a small difference exists between both methods.The numeric measures as well as the histograms were calculated on the original full size Landsat-8 scene.

Sentinel-2
The numeric assessment for the Sentinel-2 scene is in the trend with the results for the Landsat-8 data.The Sentinel-2 scene ID S2A OPER MSI L1C TL MTI 20151112T100721 A002034 T36PZS was acquired on the 12 November 2015 at 10:07:21 UTM.

CONCLUDING REMARKS
A new method for cirrus removal without the use of the 1.38µm band data has been developed.The method is based on the calculation of the cirrus thickness map and named as the CTM.The method restores the data with spectral consistency and demonstrates a small difference to the standard 1.38µm decirrus method.The developed method is fast, parameter free, and was successfully tested on a large set of Landsat-8 and Sentinel-2 data.It should be noted that an overcorrection in the 1.38µm based decirrus method can appear especially in case of moderately thick cirrus (cirrus band reflectance higher 2.5%).The presented numeric evaluation demonstrates that there is no over-/under-correction in the CTM method.In case of an absolutely intransparent cirrus none of the decirrus methods can restore the data since only a negligibly small fraction of direct radiance pass throug the intransparent cirrus cloud.The method is implemented in the AT-COR software (ATCOR, Atmospheric / Topographic Correction for Satellite Imagery.ReSe Applications Schläpfer, Langeggweg 3, CH-9500 Wil, Switzerland, n.d.) and can run in a fully automatic batch mode.Chepfer, H., Goloub, P., Shinhirne, J., Flamant, P. H., Lavorato, M., Sauvage, L., Brogniez, G., and Pelon, J., 2000.Cirrus cloud properties derived from polder-1/adeos polarized radiances: first Figures 4 b) and d) illustrate estimated cirrus map (yellow) with a slight underestimation in the 1.38µm based decirrus method.A visual analysis of the two decirrus methods (Figures 4 c), e)) reveals a slight cirrus undercorrection in Figure 4 c) (1.38µm based decirrus) compared to 4 e) (CTM decirrus).This undercorrection can be partly attributed to the underestimation of the cirrus area (compare cirrus maps in Figures 4 b) and d)).The histogram in Figure 5 a) demonstrates the mean shift for the decirrus results up to 1% and Table2demonstrates the mean shift equal to 0.56%.The peak absolute difference for the two decirrus results is up to 2.5% (Figure5 b)).The mean relative difference value is 8% (Table2, also Figure5 c)) and this value is low.Figure 6 again demonstrates a spectral consistency for the 1.38µm and CTM decirrus.As in the Landsat-8 case the 1.38µm based decirrus compensates slightly more in the shortwave infrared 1.6 and 2.2µm bands which is expected.

Table 1 :
Landsat-8 scene TOA reflectance absolute and relative difference mean value and the RMSE for each band.

Table 2 :
Sentinel-2 absolute and relative difference TOA reflectance mean value and the RMSE for each band.

Table 2
demonstrates the Sentinel-2 scene per band evaluation.Low mean absolute and relative difference values together with a low RMSE demonstrate a close agreement for the two decirrus methods.