Evaluation of Landsat-8 and Sentinel-2A aerosol optical depth retrievals across Chinese cities and implications for medium spatial resolution urban aerosol monitoring.

In urban environments, aerosol distributions may change rapidly due to building and transport infrastructure and human population density variations. The recent availability of medium resolution Landsat-8 and Sentinel-2 satellite data provide the opportunity for aerosol optical depth (AOD) estimation at higher spatial resolution than provided by other satellites. A year of 30 m Landsat-8 and 10 m Sentinel-2A AOD data retrieved using the Land Surface Reflectance Code (LaSRC) were compared with coincident ground-based Aerosol Robotic Network (AERONET) Version 3 AOD data for 20 Chinese cities. Stringent selection criteria were used to select contemporaneous data - only satellite and AERONET data acquired within 10 minutes were considered. The average satellite retrieved AOD over a 1470 m × 1470 m window centered on each AERONET site was derived to capture fine scale urban AOD variations. AERONET Level 1.5 (cloud-screened) and Level 2.0 (cloud-screened and also quality assured) data were considered. For the 20 urban AERONET sites in 2016 there were 106 (Level 1.5) and 67 (Level 2.0) Landsat-8 AERONET AOD contemporaneous data pairs, and 118 (Level 1.5) and 89 (Level 2.0) Sentinel-2A AOD data pairs. The greatest AOD values (>1.5) occurred in Beijing, suggesting that the Chinese capital was one of the most polluted cities in China in 2016. The LaSRC Landsat-8 and Sentinel-2A AOD retrievals agreed well with the AERONET AOD data (linear regression slopes > 0.96; coefficient of determination r2 > 0.90; root mean square deviation < 0.175) and demonstrate that the LaSRC is an effective and applicable medium resolution AOD retrieval algorithm over urban environments. The Sentinel-2A AOD retrievals had better accuracy than the Landsat-8 AOD retrievals, which is consistent with previously published research. The implications of the research and the potential for urban aerosol monitoring by combining the freely available Landsat-8 and Sentinel-2 satellite data are discussed.


Introduction
In much of Asia, rapid economic growth has resulted in increased air pollution [1][2][3][4]. In particular, the increased occurrence of urban haze, characterized by high concentrations of fine particulate matter with diameter less than 2.5 µm (PM2.5), has received growing attention due to human health and quality of life concerns [5,6]. In situ air quality measurements provide high temporal resolution Level 1.5 (cloud-screened) and Level 2.0 (cloud-screened and also quality assured) data were used. The AERONET AOD has a +0.02 bias and one sigma uncertainty of 0.02 [29].
All the Version 3.0 AERONET data over mainland China, Hong Kong, and Taiwan available in 2016 were considered. A whole year of data was used to capture a representative range in AOD across China. Two AERONET sites were removed however as they were not in or close to urban areas (checked by comparison with Google Map satellite images). The majority of the remaining 26 sites had both Level 1.5 and Level 2.0 data and were in mainland China (Figure 1). Site differences in the annual amount of AERONET data and the processing level (Table 1) reflect cloud conditions, and AERONET site operational differences.
Remote Sens. 2019, 11 FOR PEER REVIEW 3 higher quality Level 1.5 (cloud-screened) and Level 2.0 (cloud-screened and also quality assured) data were used. The AERONET AOD has a +0.02 bias and one sigma uncertainty of 0.02 [29]. All the Version 3.0 AERONET data over mainland China, Hong Kong, and Taiwan available in 2016 were considered. A whole year of data was used to capture a representative range in AOD across China. Two AERONET sites were removed however as they were not in or close to urban areas (checked by comparison with Google Map satellite images). The majority of the remaining 26 sites had both Level 1.5 and Level 2.0 data and were in mainland China (Figure 1). Site differences in the annual amount of AERONET data and the processing level (Table 1) reflect cloud conditions, and AERONET site operational differences.

Landsat-8 and Sentinel-2A Data
Landsat-8 was launched in 2013 and carries the optical wavelength Operational Land Imager (OLI) and the thermal wavelength Thermal Infrared Sensor (TIRS) [30]. The OLI is a 12-bit multispectral sensor providing eight reflective wavelength 30 m observations from 435 nm to 2294 nm [31]. The OLI data are provided in approximately 185 × 180 km images defined in Worldwide Reference System (WRS) path/row coordinates in the Universal Transverse Mercator (UTM) projection. Landsat-8 has Remote Sens. 2019, 11, 122 5 of 14 a 16-day revisit cycle and so nominally there are 22 or 23 OLI images over each path/row location per year. Only Landsat-8 Collection 1 Tier 1 data were used in this study as they are the highest quality data and are radiometrically calibrated and orthorectified [32]. The image-to-image registration accuracy of the Landsat-8 Collection 1 Tier 1 scenes is ≤12 m radial root mean square error [33].
Sentinel-2A was launched in June 2015 and carries the optical wavelength multi-spectral instrument (MSI). The MSI is a 12-bit multispectral sensor with four 10 m, six 20 m and three 60 m reflective bands observations from 433 nm to 2280 nm [34]. Sentinel-2A has a 10-day revisit cycle and the data are provided in fixed 109 × 109 km spatially overlapping tiles in the UTM projection [35]. The image-to-image registration accuracy of the data was 16 m before 15 June 2016 and was 4 m after this date [36] when the European Space Agency updated the processing software [37].
All of the Landsat-8 OLI and Sentinel-2A MSI data available in 2016 over the 26 urban AERONET sites were obtained. Top of atmosphere (TOA) reflectance OLI images were obtained from the United States Geological Survey portal (https://earthexplorer.usgs.gov/) and TOA Sentinel-2A L1C tile products were obtained from the European Space Agency portal (https://scihub.copernicus.eu/dhus). This provided between 14 and 66 (mean = 20) Landsat-8 OLI images and between 14 and 122 (mean = 30) Sentinel-2A tile acquisitions at each AERONET site (Table 1). More Sentinel-2A acquisitions were available because of the higher 10-day satellite revisit cycle.

Landsat-8 and Sentinel-2A Aerosol Optical Depth (AOD) Retrieval
Reliable surface monitoring with optical wavelength remotely sensed data requires atmospheric correction to minimize the scattering and absorbing effects of atmospheric gases and aerosols. Radiative transfer algorithms and atmospheric characterization data are used for automated large-area atmospheric correction [38,39]. The correction of aerosol effects is particularly challenging because aerosols are highly variable in space and time [28,40]. In this study, the LaSRC (Land Surface Reflectance Code) V3.5.5 atmospheric correction algorithm [41] was used. LaSRC code is publicly available (https://github.com/USGS-EROS/espa-surface-reflectance/tree/master/lasrc) and has been operationally used by the United States Geological Survey (USGS) to generate Landsat Analysis Ready Data (ARD) [32] and by NASA to generate Harmonized Landsat and Sentinel-2 surface reflectance products [42]. The LaSRC algorithm was developed for atmospheric correction of Landsat-8 imagery [41] and has been adapted for Sentinel-2A application [26,43]. It is based on the 6SV radiative transfer code [44]. The AOD retrieval takes advantage of the short wavelength red and blue bands on the Landsat-8 and Sentinel-2 sensors. Two surface reflectance ratios, red to blue and red to ultra-blue, and the difference between them are used to invert the AOD using a global coarse resolution (0.05 • ) ratio data set derived from MODIS and MISR data, and expressed as a function of a mid-infrared vegetation index [41]. A fixed urban-clean aerosol model [28] is assumed. Despite using a fixed aerosol model, the Ångström exponent (related to the dependence of AOD with wavelength) is not fixed [41]. This has a similar effect as allowing a dynamic aerosol model with varying aerosol size distribution and refractive index [14,40] and is needed for aerosol retrieval over regions with different aerosol types. The LaSRC AOD is defined at 550 nm for each 10 m (Sentinel-2A) and 30 m (Landsat-8) pixel, and the LaSRC code was modified to write out the AOD values. The LaSRC algorithm also generates a per-pixel cloud mask using the red-to-blue surface reflectance ratios [41]. The LaSRC 10 m (Sentinel-2A) and 30 m (Landsat-8) cloud masks were used in this study to remove cloud contaminated pixels.

Comparison of Contemporaneous AERONET and Satellite Retrieved AOD
The ground-based AERONET 500 nm AOD measurements were compared with the LaSRC satellite retrieved AOD 550 nm data. First, the AERONET 500 nm AOD measurements were converted to 550 nm equivalent values using the standard Ångström exponent interpolation method [45] as: where τ(550) is the interpolated AERONET 550nm AOD data, τ(500) is the AERONET 500 nm AOD measurement, and α is the AERONET 440-675 nm Ångström exponent. Other studies have suggested that a quadratic or cubic relationship may better characterize the AOD dependence with wavelength in the logarithmic scale [46]. However, as the interpolated wavelength (550 nm) is close to the observed wavelength (500 nm), the log-linear relationship is used as Equation (1), and following on from well-established MODIS approaches [47,48].
The following temporal and spatial data selection criteria were applied to select contemporaneous satellite and AERONET AOD data. At each AERONET site the closest AERONET AOD 550 nm value within 10 min (before or after) of the satellite overpass time was selected, and the mean 550 nm LaSRC AOD over a 1470 m × 1470 m image window centered on the AERONET site location was derived. If more than half of the satellite AOD retrievals in the 1470 m × 1470 m window were cloudy or missing, the data were discarded. A 10-min period was used because the median closest temporal difference between the satellite overpass times and the AERONET measurements for the 2016 data was 2.9 min (Sentinel-2A) and 4.1 min (Landsat-8), and the median second closest temporal difference was 7.7 min (Sentinel-2A) and 9.5 min (Landsat-8). In China, the annual urban mean wind speed is 2.4 m/s [49] and so in 10 min aerosols could be blown 1440 m, which is smaller than the 1470 m window dimension. Using too large a window size will increase the likelihood of averaging aerosols from different sources, which is a concern in urban environments where aerosols may vary spatially quite rapidly. Table 2 summarizes the number of contemporaneous satellite and AERONET AOD data sets in 2016. Due to the stringent selection criteria (Section 3.2), six of the 26 AERONET sites had no remaining contemporaneous data, leaving a total of 20 urban AERONET sites considered in the reminder of this study. There were generally more contemporaneous Sentinel-2A than Landsat-8 data due to the greater Sentinel-2A temporal data availability. There were more Level 1.5 data than Level 2.0 AERONET data due to the additional quality assurance used to generate the Level 2.0 data [29]. On average, among the 20 sites in 2016, there were 5.9 Level 1.5 and 4.45 Level 2.0 AERONET AOD and Sentinel-2A contemporaneous data pairs, and on average for Landsat-8, there were 5.3 Level 1.5 and 3.35 Level 2.0 AERONET AOD data pairs.

Contemporaneous Data Availability and Example Annual Satellite and AERONET AOD Data Comparison
The Beijing-CAMS and XiangHe AERONET sites had the greatest number of contemporaneous satellite and AERONET data pairs. AOD retrieval [50,51]. Examination of the other days with contemporaneous data in 2016 indicated no other cloud detection omission errors. In summary, for the illustrated year of Beijing-CAMS site data, after discarding the day 181 cloud contaminated data, the root mean square deviation (RMSD) between the AERONET AOD and the satellite AOD was 0.183 for Landsat-8 (n = 12) and 0.102 for Sentinel-2A (n = 28). The correlation between the satellite AOD and the AERONET AOD for the year was 0.967 for Landsat-8 and 0.962 for Sentinel-2A. Table 2. Summary of the 26 urban AERONET sites (Figure 1), the number of days with 2016 AERONET data, and the number of contemporaneous LaSRC satellite AOD and AERONET AOD data pairs.  AOD retrieval [50,51]. Examination of the other days with contemporaneous data in 2016 indicated no other cloud detection omission errors. In summary, for the illustrated year of Beijing-CAMS site data, after discarding the day 181 cloud contaminated data, the root mean square deviation (RMSD) between the AERONET AOD and the satellite AOD was 0.183 for Landsat-8 (n = 12) and 0.102 for Sentinel-2A (n = 28). The correlation between the satellite AOD and the AERONET AOD for the year was 0.967 for Landsat-8 and 0.962 for Sentinel-2A.     Notably, more than 13% of the AERONET Level 2.0 AOD values were >1.0 (indicating very hazy urban conditions). For five overpass dates the AERONET AOD values were >1.5; all occurred in Beijing (at AERONET sites: Beijing, Beijing-CAMS, and XiangHe) in the months of March, June, July, August, and September. These high AOD values were not due to undetected clouds, as the AERONET Ångström exponent values for the five overpass dates varied from 0.8680 to 1.3127 (mean 1.0941) and clouds typically have Ångström exponent values <0.5 [52]. The greatest AERONET AOD value with contemporaneous non-cloudy Landsat-8 satellite data was 1.940 and occurred at the Beijing AERONET site on July 7th (day 189), and the greatest AERONET AOD value with contemporaneous non-cloudy Satellite-2A satellite data was 2.094 and occurred at the Beijing_PKU site on August 11th (day 224). On these two days the difference between the satellite and AERONET AOD values was <0.12. These results, and those illustrated in Figures 3 and 4, indicate the ability of the LaSRC satellite AOD retrieval algorithm to work over a wide range of urban AOD values.

Discussion
The criteria used to select contemporaneous medium resolution satellite and AERONET data across China were purposefully stringent to capture fine scale urban aerosol variations. The closest AERONET AOD value within 10 min of each satellite overpass was selected, and the mean LaSRC AOD over a 1470 m × 1470 m image window centered on the AERONET site location was derived. A ±10 min period was used to reflect the approximately three minutes frequency of most AERONET data [29] and because the median closest and second closest temporal differences between the satellite overpass times and the AERONET measurements for the 2016 data were less than 10 min. A 1470 m × 1470 m image window size was used, as it corresponds to an integer multiple of the Landsat-8 30 m and the Sentinel-2A 10 m bands, and because in 10 min aerosols would not be transported completely across the window at the reported annual urban 2.4 m/s mean wind speed [49]. We undertook a sensitivity analyses, considering 2970 m × 2970 m and 5970 m × 5970 m window dimensions for the same ±10 min period, and found only small differences in the reported results (differences only in the second decimal place of the reported RMSD values).
Using large window sizes and long temporal periods as selection criteria will increase the likelihood of aerosols from different sources being compared between the AERONET and satellite AOD data. This is a concern in urban environments where aerosols may vary spatially quite rapidly. For example, Chen et al. [53] considered a year of PM2.5, PM10, and greenhouse gas concentration data from 35 Beijing municipal environmental monitoring stations and found significant spatial differences between urban, suburban, and traffic sites, and significant diurnal variations at individual sites. Similarly, Wang et al. [54] documented urban PM2.5 concentrations in the morning and evening that were nearly twice those in the mid-afternoon in northeast China. Previously, researchers have compared coarser spatial resolution but near daily MODIS AOD retrievals with AERONET AOD using temporal selection criteria of ±30 min and spatial window side dimensions varying from 9 km to 50 km [9,20,55-58] and ±7.5 min with a 10 km spatial window dimension [59]. An explicit study of spatial scaling effects, comparing MODIS AOD data with PM2.5 measurements in the Boston metropolitan area [22], found reduced correlation with increasing window size. Further research to undertake scaling analyses, including wind speed information, using the Landsat-8 and Sentinel-2 AOD data is recommended, subject to the availability of data from a dense ground-based atmospheric monitoring network.
In this study ordinary least squares (OLS) linear regression was used to fit the LaSRC satellite AOD against the AERONET AOD. This is the conventional approach adopted by the majority of researchers because the errors in the AERONET AOD data are considerably smaller than those typically present in satellite AOD retrievals. We also examined the use of reduced major axis regression (RMA) which allows error in both the independent and dependent variables [60,61]. However, the OLS and RMA results were similar with linear regression line slope differences <0.05.
Recently, the LaSRC AOD algorithm was assessed, with seven other AOD algorithms, considering a year of Landsat-8 and seven months of Sentinel-2A overpasses at 19 globally distributed urban and non-urban AERONET sites using a ±15 min and 9 km window dimension [26]. The authors [26] reported that the Sentinel-2A LaSRC AOD performed the best among the different AOD algorithms, and, as in this study, the Landsat-8 LaSRC AOD was less accurate than the Sentinel-2A LaSRC AOD. The reasons for this sensor difference are unknown but are likely related to the higher spatial resolution of the Sentinel-2A compared to Landsat-8 and to sensor spectral band differences. The uncertainty of the LaSRC derived AOD is due to a number of sources, like other AOD retrieval algorithms ( [8,41,62]), including the sensor radiance calibration uncertainty, cloud detection omission errors, and assumptions concerning the aerosol type and spectral variation in surface reflectance.
The current state-of-the-practice for global coverage satellite AOD monitoring is based on near daily, but coarse spatial resolution, polar-orbiting data, such as from MODIS, that may be spatially too coarse for monitoring urban pollution [22,63,64]. The Landsat-8 and Sentinel-2A satellites provide a higher spatial resolution (10 m to 30 m) global coverage AOD monitoring capability but with lower temporal resolution. For example, in this study there were an average of 20 Landsat-8 and 30 Sentinel-2A overpasses per year over the Chinese cities. Data from the recently launched (March 2017) Sentinel-2B were not available for the study period, but combined with Sentinel-2A and Landsat-8, will increase the temporal observation frequency to a global median average revisit interval of 2.9 days with different overpass times [65]. Given their high accuracy in AOD retrieval, future research to investigate the combined use of Landsat-8 and the two Sentinel-2 satellites for AOD monitoring is recommended.

Conclusions
In 2016 there were 26 urban AERONET sites across China with available data, and after applying the spatial and temporal selection criteria 20 sites remained with an annual total of 106 (Level 1.5) and 67 (Level 2.0) Landsat-8 AERONET AOD contemporaneous data pairs, and 118 (Level 1.5) and 89 (Level 2.0) Sentinel-2A AOD data pairs. The resulting AERONET data encompassed a wide range of AOD from clear days (minimum = 0.008) to very hazy conditions (maximum = 2.094). The greatest AOD values (>1.5) occurred in Beijing, suggesting that the Chinese capital was one of the most polluted cities in China in 2016.
The Landsat-8 and Sentinel-2A AOD retrievals agreed well with the AERONET AOD measurements (OLS linear regression slopes > 0.96, r 2 > 0.90, and RMSD < 0.175) with no significant pattern of over-or under-estimation of the AOD. Although the AERONET Level 2.0 AOD data are quality assured, there were no pronounced differences between the results considering the Level 1.5 and considering the Level 2.0 AOD data. The Sentinel-2A AOD retrievals had slightly better agreement (RMSD < 0.11) with the contemporaneous AERONET AOD measurements than the Landsat-8 retrievals (RMSD < 0.175).