Aerosol optical depth (AOD) retrieval for atmospheric correction in Landsat-8 imagery using second simulation of a satellite signal in the solar spectrum-vector (6SV)

Atmospheric correction has been challenging task in digital image processing. It requires several atmospheric parameters in order to obtain accurate surface reflectance of objects within the image scene. One of the most crucial parameters required for accurate atmospheric correction is aerosol optical depth (AOD). AOD can be obtained by in-situ measurement or estimated from remote sensing observation. In this experiment, atmospheric correction was performed using second simulation of a satellite signal in the solar spectrum-vector (6SV) algorithm on Landsat-8 imagery in which AOD parameter was retrieved from surface reflectance inversion involving daily-global surface reflectance product of moderate resolution imaging spectroradiometer (MODIS). Furthermore, AOD retrieved from surface reflectance inversion was also validated using ground-based sun photometer observation data from aerosol robotic network (AERONET) station in Bandung, Indonesia. Our experiment shows the consistency between AOD from surface reflectance inversion and AOD from ground-based observation. Finally, 6SV was performed on Landsat-8 imagery to obtain the surface reflectance. We further compared surface reflectance of 6SV atmospheric correction and surface reflectance of Landsat-8 Level 2 product. The atmospherically corrected image also shared agreeable result with Landsat 8 Level-2 product.


Introduction
Atmospheric correction has been a real challenge in preimage processing. Atmospheric effects due to aerosol in particular can be found extremely difficult to correct [1]. Aerosols are basically liquid or solid particles suspended in the atmosphere with variety of sizes, shapes and compositions [2]. Aerosols vary temporally and spatially. It affects Earth's climate systems as it plays an essential role in Earth's solar radiation budget [3]. In remote sensing, the presence of aerosols in the atmosphere ultimately affect the recorded data by spacebased instruments due to its solar radiation scattering and absorption properties.
A number of measurements and experiments have been conducted to study and measure aerosols optical depth (AOD) in the atmosphere [4]- [6]. AOD is the measure of scattering and absorption of solar radiation due to aerosol particles [3]. AOD is one of the compendious variable to assess using ground-based instrument which represents the characteristics of local aerosol conditions such as aerosol pollution [7]. AOD is also used to make atmospheric correction in remote sensing data. However, ground-based aerosol observation data are still limited due to the limitation of the instruments. Aerosol optical depth retrieval using space-based instruments such as from Moderate Resolution Imaging Spectroradiometer (MODIS) product provides opportunity to estimate AOD parameter for atmospheric correction application.
Several methods of atmospheric correction have been proposed and implemented such as 6SV, FLAASH, QUAC, ATCOR, etc. These methods can be applied to remove atmospheric effects in many different optical satellite imageries, including Landsat imagery. The latest satellite of Landsat mission, Landsat-8, was launched in 2013. Landsat-8 data have been widely used in various fields and applications such as land cover monitoring and bathymetry mapping [8], [9].
In some applications, Landsat-8 surface reflectance (SR) is prerequisite for further analysis such as in land cover monitoring and bathymetry extraction [10], [11]. The surface reflectance product is the reflectance value of object on the land without atmospheric effects [12] . Atmospheric effects are eliminated by performing atmospheric correction in order to obtain SR. Standardised Landsat data processing for SR production has been developed and also demonstrated by Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) project and Web-enabled Landsat Data (WELD) project [13], [14]. The overall objective of both projects is to characterise and monitor land cover changes using SR product of Landsat imageries. The main difference between LEDAPS and WELD is aerosol characterisation as input to perform atmospheric correction. Aerosol parameter is derived independently from each Landsat acquisition using ancillary information of water vapour, ozone and barometric pressure in LEDAPS processing. However, WELD uses Terra-MODIS products to correct atmospheric effects in Landsat Enhanced Thematic Mapper (ETM+) imagery [15].
The objectives of this study are (1) to demonstrate AOD retrieval using MODIS products and (2) to perform atmospheric correction using 6SV in Landsat-8 imagery based on the retrieved AOD parameter. We also compared SR of the corrected image and Landsat-8 Level 2.

Landsat-8 Data
In this experiment, Landsat-8 Level 1 product (path/row: 122/65) which covers Kota Bandung, Indonesia was used to be atmospherically corrected (Fig.1a). Landsat-8 Level 1 product is represented by digital number values which can be converted to Top-of-Atmosphere reflectance (TOA). Moreover, Landsat-8 Level 2 SR Product was also used to compare SR value obtained from atmospheric correction. Landsat-8 Level 2 ( Fig.1.b) SR was used as reference to evaluate atmospheric correction based on AOD inversion in this experiment. Both datasets were recorded on 28 June 2015. While Landsat-8 Level 1 product can be downloaded directly from https://earthexplorer.usgs.gov/ , Landsat-8 Level 2 SR product can be acquired on demand at no cost from the mentioned website.

MODIS surface reflectance product
MODIS has been providing a daily and nearly global scale remote sensing data using Terra and Aqua instrument on board for more than 10 years. MODIS allows scientist to observe entire planet to have better understanding by providing remote sensing data such as atmospheric parameters and including SR product [16]. In this experiment, we utilised Aqua-MODIS SR product with coarse spatial resolution of 0.05° (Fig.2) to perform AOD inversion along with Landsat-8 Level 1 product. Aqua-MODIS surface reflectance product was acquired on the same day with Landsat-8 imagery used. Aqua-MODIS SR product was used as primary data to retrieve AOD parameter in Landsat-8 alternative to Terra-MODIS SR product [1], [17], [18]. This alternative was taken because of the gapped region right over the study area in Terra-MODIS product. Therefore, Aqua-MODIS SR product of MDY09CMG and MYD09CMA are used. SR products of band 1, 2, 3 and 7, and 8 are used to estimate AOD parameter. The characteristic and specification of each band is presented in Table 1. SR data of band 1, 2, 3 and 7 are found in MDY09CMG product, while Aqua SR of band 8 is found in MYD09CMA product which has to be extracted in order to be processed with other bands for AOD inversion. MODIS bands specification used in this study is presented in Table 1.

AOD ground-based observation data
AOD ground-based observation data were acquired from Aerosol Robotic Network (AERONET). AERONET is a global network of distributed ground-based stations that provide daily aerosol optical thickness or AOD data. (Ju et al., 2012) There are 8 AERONET stations in Indonesia, however Bandung AERONET station is the one that provides near-continuous AOD data. We used AOD Level-2 based on Version 2-Direct Sun Algorithm which was acquired from AERONET website. The data included AOD observations from 12 May 2009 to 13 February 2017. However, the instrument only measured and observed 10 data throughout the day on 28 June 2015.

Landsat-8 Radiometric Calibration
Landsat-8 Level 1 product is presented by 16-bit digital number (DN) on each pixel. The digital number can be converted to radiance (L λ ) or TOA reflectance (ρλ). Radiometric calibration is performed to obtain L λ or ρλ with the provided parameters in Landsat-8 Level. Radiance (L λ ) or TOA reflectance (ρλ) can be calculated as follows: where M L is spectral radiance multiplicative scaling factor, A L is radiance additive scaling factor, M ρ is reflectance multiplicative scaling factor, A ρ is reflectance additive scaling factor, θ is solar angle and Q cal is DN.

AOD Inversion
Coastal/aerosol band is one of Landsat 8 OLI channels which has wavelength range of 433 -453 nm. This channel has two main functions, namely 1) to estimate aerosol concentrations in the atmosphere to be used for atmospheric corrections, and 2) to be used in to coastal zone studies. AOD describes the influence of aerosols to the transmission of light. In this experiment, AOD retrieval method is primarily adapted from previous study [1] in which AOD inversion is performed to retrieve AOD in Landsat-8 using surface reflectance product of Aqua-MODIS. Coastal/aerosol band of OLI was also used to invert the AOD. The basic principle of AOD inversion relies on the relationship between the reflectance in MODIS middle infrared (Mid-IR) band (2.1 nm) and coastal/aerosol band (0.43 nm). The relationship is based on the fact that liquid water content from photosynthesis of vegetation absorbs solar radiation in coastal/aerosol band and also affects the observed reflectance on Mid-IR band. This inversion utilises Normalised Difference Vegetation Index (NDVI) using NIR band (band 2) and Mid-IR (band 7) instead of red band in Aqua-MODIS. NDVI MIR is computed with the following equation: It is worth mentioning that in this experiment, the presence of thin cloud or cirrus cloud was not taken into account, neither water vapour, ozone nor barometric pressure were considered.

AOD Validation using AERONET Data
The AOD-retrieved from OLI is validated using AERONET observation data on the same day when OLI imagery was recorded. AOD-retrieved from Aqua-MODIS is measured at wavelength of 470, 550 and 660 nm, while AEORNET data is derived at wavelength of 340, 380, 440, 500, 675, 870 and 1020 nm in available data. Moreover, due to the limitation of AERONET observation data on 28 June 2015, it is not possible to do band interpolation for particular band.
We recognise unideal conditions within this experiment due to limited AERONET data. However, comparison was still made between AOD-retrieved at 550 nm from OLI and MODIS with AOT AERONET at 500 nm. We calculated expected error (EE) based on the EE equation on previous study [19]. EE is calculated to evaluate the quality of AOD-retrieved in this experiment. EE can be expressed as follows: EE = ± (0.05 + 0.15AOD AERONET ) (4) AOD-retrieved is considered to have good quality when AODretrieved falls within EE envelope: AOD AERONET − |EE| ≤ AOD ≤ AOD AERONET + |EE| (5)

6SV Atmospheric Correction
The 6S is radiative transfer code used for the calculation of look-up tables in MODIS atmospheric correction. It is among the most commonly used, strictly validated with comprehensive documented radiative transfer code in remotesensing community. Atmospheric correction process is carried out to retrieve surface reflectance based on simulations of atmospheric effects such as molecular scattering, gaseous absorption and aerosol that ultimately affect TOA reflectance. TOA reflectance can be expressed as follows [12]: where ρ TOA is TOA reflectance, T g is to gaseous transmission, ρ R is the molecular scattering intrinsic reflectance, ρ R+A is scattering intrinsic reflectance of molecules and aerosols, T R+A refers to transmission of molecules and aerosols, S R+A is spherical albedo and ρ s is surface reflectance.
Atmospheric correction is made by using internet-based 6SV software on http://6s.ltdri.org/. Surface reflectance can be retrieved based on atmospheric correction coefficients of X A , X B and X C . SR is expressed by following equation [20]: where is SR value, ( X A , X B , X C ) is 6SV correction coefficients and L λ refers to OLI radiance value. 6SV requires several parameters as inputs to obtain atmospheric correction coefficients. Atmospheric correction coefficients were used to eliminate atmospheric effect from TOA reflectance to obtain SR. In this experiment, we used several options as input parameters which are presented in Table 2.

AOD Retrieval
AOD retrieval in Landsat-8 Level 1 imagery was carried out by utilising Aqua-MODIS SR product. In this experiment, we have successfully retrieved AOD (550nm) for entire OLI scene (path/row: 122/65) which ranges from 0.06978 to 0.911947. Furthermore, image cropping was made based on administrative boundary of Kota Bandung to obtain AOD value in Kota Bandung. It is found that AOD in Kota Bandung of OLI ranges from 0.118213 to 0.911947 (Fig.3)

Fig. 3. AOD map over Kota Bandung
AOD map shows a pattern which follows land cover over study area. Based on the AOD map, vegetation area shows lower AOD than its surrounding. On the other hand, it is found that roads, urban and industry area tend to have higher AOD. Previous studies suggest that human activities contributes to high AOD [21], [22].

AOD Validation
As mentioned earlier, there are several unideal conditions in this experiment and mainly due to the limited ground-based AOD measurements. Ideally, well-distributed AERONET stations with near-continuous AOD data are required in order to accurately validate satellite-retrieved AOD parameters. We also recognise that it is not reasonable to compare single pixel AOD-retrieved value to AERONET measurement point. However, the available AERONET data only provides one possible station to be used and even that, the limitation of observation data is still an obstacle in this experiment.
AERONET station in Bandung only measured and observed 10 data of AOD on 28 June 2015 which included 8 data observed around 5 hours before Landsat-8 passing time and 2 data observed at midnight. Landsat-8 passing time over Bandung is at 09:59:50 AM, while AERONET AOD data which is the closest to Landsat-8 passing time was observed at 04:32:47 AM. Therefore, there was huge gap between Landsat-8 passing time and AERONET observation time.
Based on Equation (4) and (5), we evaluated the quality of single pixel AOD-retrieved AOD with AERONET data. We found that AOD-retrieved in this experiment falls within the EE. AOD-retrieved of 0.209521 falls within the range of 0.193023 ≤ AOD-retrieved ≤ 0.378795. Although AODretrieved can be considered as good statistically, again it is not realistic to compare AOD observation data at 04:32:47 AM to AOD-retrieved value at 09:59:50 AM since the atmosphere is constantly changing and aerosols are varying temporally.

SR Comparison
Atmospheric correction was carried out to eliminate atmospheric effect in TOA reflectance and to obtain surface reflectance value of Landsat-8 imagery using 6SV algorithm. In this experiment, we did not perform atmospheric correction on the entire scene to obtain SR value for each pixel, rather we only performed atmospheric correction for 15 pixels as sample which included 5 pixels for each vegetation, urban and waterbody. Atmospheric correction was carried out in band 1 and band 4 OLI to evaluate the sensitivity of both spectral bands used in AOD inversion. The comparison between SRretrieved from 6SV and SR-product from Landsat-8 Level 2 for band 1 is presented in Figure 4.
Based on the comparison and difference (Δ SR ) between two SR products, Root Mean Square Error (RMSE) was calculated for all of 15 pixels and resulted in RMSE of 0.015. The comparison suggests that water body pixels do not show strong agreement between two SR products. Correlation between two SR-products is also made and presented as linear function (Fig.  5). It is shown that the linear regression between two SR products does not have good quality which is represented by the value of R 2 . The R 2 of 0.404 is considered low which means that both SR datasets do not have a strong correlation. However, it is worth mentioning that SR comparison in coastal aerosol band is not idealistic to make due to the fact that this band was used as input to retrieve AOD parameter. 6SV atmospheric correction in Landsat-8 band 4 also retrieved SR-product which was compared to Landsat-8 Level 2 SR product (Fig.4). Atmospheric correction for OLI band 4 gives RMSE of 0.013 which means that the results have better quality than 6SV done in OLI band 1. Correlation between both SR products in OLI band 4 also shows consistency with the R 2 of 0.9197 (Fig.6). Therefore, we consider that the atmospheric correction in OLI band 4 was performed quite well with strong agreement between SR-retrieved and Landsat-8 Level 1 product as reference.

Conclusion
This study has demonstrated AOD inversion using Aqua-MODIS product to retrieve AOD map on Landsat-8 OLI imagery. OLI AOD-retrieved is further used as an input to atmospherically correct OLI imagery to obtain surface reflectance product using 6SV atmospheric correction. There are several limitations in this experiment which are the limitation of ground-based instrument to measure local AOD and limited available AOD data. However, with many unideal conditions, our experiment has successfully retrieved AOD map for Landsat-8 which falls within the expected error. Furthermore, 6SV correction also shows high consistency between SR-retrieved with SR-product generated from Landsat-8 Level 2. It is worth mentioning that our approach is not suitable for monitoring applications such as AOD retrieval for air quality monitoring and SR production for land cover monitoring. This experiment is merely a practical demonstration to retrieve AOD in Landsat-8 and to perform atmospheric correction based on AOD-retrieved.