Integration of Surface Reflectance and Aerosol Retrieval Algorithms for Multi-Resolution Aerosol Optical Depth Retrievals over Urban Areas

The SEMARA approach, an integration of the Simplified and Robust Surface Reflectance Estimation (SREM) and Simplified Aerosol Retrieval Algorithm (SARA) methods, was used to retrieve aerosol optical depth (AOD) at 550 nm from a Landsat 8 Operational Land Imager (OLI) at 30 m spatial resolution, a Terra-Moderate Resolution Imaging Spectroradiometer (MODIS) at 500 m resolution, and a Visible Infrared Imaging Radiometer Suite (VIIRS) at 750 m resolution over bright urban surfaces in Beijing. The SEMARA approach coupled (1) the SREM method that is used to estimate the surface reflectance, which does not require information about water vapor, ozone, and aerosol, and (2) the SARA algorithm, which uses the surface reflectance estimated by SREM and AOD measurements obtained from the Aerosol Robotic NETwork (AERONET) site (or other high-quality AOD) as the input to estimate AOD without prior information on the aerosol optical and microphysical properties usually obtained from a look-up table constructed from long-term AERONET data. In the present study, AOD measurements were obtained from the Beijing AERONET site. The SEMARA AOD retrievals were validated against AOD measurements obtained from two other AERONET sites located at urban locations in Beijing, i.e., Beijing_RADI and Beijing_CAMS, over bright surfaces. The accuracy and uncertainties/errors in the AOD retrievals were assessed using Pearson’s correlation coefficient (r), root mean squared error (RMSE), relative mean bias (RMB), and expected error (EE = ± 0.05 ± 20%). EE is the envelope encompassing both absolute and relative errors and contains 68% (±1σ) of the good quality retrievals based on global validation. Here, the EE of the MODIS Dark Target algorithm at 3 km resolution is used to report the good quality SEMARA AOD retrievals. The validation results show that AOD from SEMARA correlates well with AERONET AOD measurements with high correlation coefficients (r) of 0.988, 0.980, and 0.981; small RMSE of 0.08, 0.09, and 0.08; and small RMB of 4.33%, 1.28%, and -0.54%. High percentages of retrievals, i.e., 85.71%, 91.53%, and 90.16%, were within the EE for Landsat 8 OLI, MODIS, and VIIRS, respectively. The results suggest that the SEMARA approach is capable of retrieving AOD over urban areas with high accuracy and small errors using high to medium spatial resolution satellite remote sensing data. This approach can be used for aerosol monitoring over bright urban surfaces such as in Beijing, which is frequently affected by severe dust storms and haze pollution, to evaluate their effects on public health.


Introduction
Aerosol optical depth (AOD) is retrieved from satellite observations using a variety of methods depending, among others, on the sensor characteristics and the retrieval approach chosen [1,2]. Satellite-retrieved AOD is used by a wide scientific community to monitor and investigate atmospheric aerosol properties at regional to global scales. AOD enables understanding of how aerosols affect the Earth's climate system [3], atmospheric visibility [4], and public health [5]. Ground-based Sunphotometer networks, such as the AErosol RObotic NETwork (AERONET) [6], have been established worldwide for continuous measurements of aerosol optical properties, including AOD, with high temporal resolution. However, these point measurements have sparse coverage. This spatial limitation is overcome by satellite remote sensing, which provides near-real-time AOD information at a global scale. Satellite remote sensing provides operational AOD products at spatial resolutions of 1 km, 3 km, 4.4 km, 6 km, and 10 km from several sensors including Moderate Resolution Imaging Spectroradiometer (MODIS), Multi-angle Imaging Spectroradiometer (MISR), Visible Infrared Imaging Radiometer Suite (VIIRS), etc., [7][8][9][10][11][12][13][14][15][16]. AOD inversion is mainly sensitive to an accurate cloud mask, an appropriate selection of the aerosol model, an accurate estimation of the surface reflectance, or an effective decoupling of the surface and atmospheric contributions to the reflectance at the top-of-atmosphere (TOA) [1,17,18]. Similarly, an accurate surface reflectance estimation based on a Radiative Transfer Model (RTM) depends on precise and accurate AOD retrievals. In other words, the retrieval of either AOD or surface reflectance requires accurate knowledge of the effects of the other parameters on the measured TOA reflectance, and thus, retrieval algorithms for these variables depend on each other.
In the present study, surface reflectance was estimated using the SREM (Simplified and Robust Surface Reflectance Estimation Method) atmospheric correction method [19,20]. SREM is the simplest atmospheric correction method compared to other existing imagebased and physical atmospheric correction methods. It performs surface reflectance inversion for VIS-SWIR bands based on the Satellite Signal in the Solar Spectrum-Vector (6SV) [21][22][23] RTM equations 'without' incorporating information on AOD, aerosol model, water vapor concentrations, ozone concentrations, and atmospheric gases. It only requires TOA reflectance/radiance and solar-sensor geometries to estimate surface reflectance. Previous studies [19,20] showed that the SREM estimated surface reflectance is similar to the surface reflectance obtained from operational surface reflectance products of Landsat 8 OLI (Operational Land Imager), MODIS (Moderate Resolution Imaging Spectroradiometer), and VIIRS (Visible Infrared Imaging Radiometer Suite). It should be noted that the SREM estimated surface reflectance is different from the Rayleigh-corrected TOA reflectance, which is only a subtraction of the Rayleigh reflectance from the TOA reflectance and does not incorporate total transmission and the atmospheric backscattering ratio. Details of the SREM atmospheric correction method can be found in Bilal et al. [19].
Monitoring of atmospheric aerosol properties over complex and mixed bright urban surfaces, such as in Beijing, is challenging, especially during high aerosol loading. The Simplified high-resolution Aerosol Retrieval Algorithm (SARA) has been developed using MODIS data with a spatial resolution of 500 m to retrieve AOD at a regional scale. SARA does not create a comprehensive look-up table (LUT) for AOD inversion and can retrieve AOD directly for the green channel. SARA requires TOA reflectance/radiance, solarsensor zenith and azimuth angles, and aerosol information to perform aerosol inversion by iterating a wide range of aerosol types (ω o = 0.30-1.0 and g = 0.0-1.0), as SARA does not have prior information of aerosol types over the region. Aerosol information can be obtained from the local AERONET site or any reliable remote sensing aerosol product. SARA results have been validated over the complex and bright urban surfaces in Hong Kong [24], over the mixed urban surfaces in Beijing during dust storms [25] and haze episodes [26], and over water surfaces [11]. The SARA retrieved AOD has also been used for the estimation of concentrations of fine particulate matter (PM 2.5 ) at 500 m resolution over Hong Kong and the Pearl River Delta (PRD) [27]. The SARA results are very promising compared to the previously published high-resolution aerosol retrieval algorithms [28][29][30].
The objective of this study is to introduce a new regional approach 'SEMARA', i.e., an integration of the SREM and SARA methods, which can accurately retrieve regional AOD over bright urban surfaces, such as Beijing, from multi-resolution remote sensing data. In previous studies of SARA AOD, surface reflectance was obtained from the operational surface reflectance products. For example, in the case of MODIS and Landsat 8 OLI, surface reflectance was obtained from MOD09 [24,25] and LaSRC [31] surface reflectance products, respectively. These operational products show positive and negative biases in surface reflectance over bright urban and hilly terrain surfaces, respectively, compared to the surface reflectance estimated by [19,20]. Therefore, the SEMARA integrated approach is expected to be an improvement as it obtains corrected surface reflectance from the SREM compared to MOD09 and LaSRC.

Satellite Data
Remote sensing images from three satellite sensors, including Landsat 8 OLI, MODIS, and VIIRS, were obtained to retrieve AOD using the SEMARA approach described in

AERONET Data
AERONET Version 3 (V3) Level 1.5 (L1.5) AOD measurements [32] were obtained for the Beijing_RADI site from 2013 to 2021. V3 L2.0 AOD data at the Beijing site from 2013 to 2017 and V3 L1.5 AOD data at the Beijing_CAMS site from 2017 to 2021 were obtained from the AERONET website (https://aeronet.gsfc.nasa.gov/; accessed on 14 August 2021) ( Table 1). In our previous study [26], we found a Pearson's correlation coefficient (r) = 1, slope = 1, and intercept = 0 between L1.5 and L2.0 AOD measurements. Therefore, L1.5 AOD was used when L2.0 AOD data were not available. AOD measurements from the Beijing site were used as an input in the SEMARA approach as explained in the research methodology section (Section 3). AOD data from the Beijing_RADI and Beijing_CAMS sites were used to validate the SEMARA retrieved AOD.

Description of the SEMARA Approach
In the present study, AOD was retrieved with spatial resolutions of 30 m, 500 m, and 750 m using data in the green channels of Landsat 8 OLI, Terra-MODIS, and SNPP-VIIRS by application of the SEMARA approach (Equation (1a,b)), an integration of the SREM (Equation (1a)) [19] and the SARA (Equation (1b)) [24] algorithms.
where ρ TOA = the top-of-atmosphere (TOA) reflectance, which is a function of measured spectral radiance (L TOA ), solar zenith angle, earth-sun distance (d) in the astronomical unit, and mean solar exoatmospheric radiation (E SUN ). ρ R = Rayleigh reflectance in the absence of aerosols. λ = wavelength. T s = atmospheric transmittance on the sun-surface path (downward). T s0 = same as T s but in an aerosol-free atmosphere.
T v = atmospheric transmittance on the surface-sensor path (upward). T v0 = same as T v but in an aerosol-free atmosphere. ρ s = SREM estimated surface reflectance. S atm = atmospheric backscattering ratio to account for multiple reflections between the surface and atmosphere. S atm = same as S atm but in an aerosol-free atmosphere. µ s = cosine of the solar zenith angle. µ v = cosine of the sensor zenith angle. P a = aerosol phase function. ω o = single scattering albedo. τ a = SARA AOD.
where λ = wavelength in µm τ R = Rayleigh optical depth θ s = solar zenith angle θ v = sensor zenith angle ϕ = relative azimuth angle A and B are coefficients that account for the molecular asymmetry.
where g is the asymmetry factor.
In Equations (5a-c) and (7), τ A (AOD, where A represents AERONET) is an unknown parameter that can be obtained from the AERONET Sunphotometer or other remote sensing aerosol products. In the present study, τ A is obtained from the Beijing AERONET site, another urban site besides Beijing_RADI and Beijing_CAMS, to solve Equations (5a-c) and (7). In Equations (1a,b), (5a-c)- (7), ω o and g are unknown. The ω o and g are derived using AERONET AOD (τ A ), which is required on the right-hand side of Equation (1b) for the calculation of the total atmospheric transmittances (T s T v ) and atmospheric backscattering ratio (S atm ). This is achieved by varying ω o and g within the ranges of 0.30-1.0 and 0.0-1.0, respectively, until AOD (τ a ) on the left-hand side (τ a ) of Equation (1b) becomes equal to the input AERONET AOD (τ A ) on the right-hand side of Equation (1b). In this way, the values of ω o and g are selected such that τ a equals τ A in Equation (1b) when the TOA reflectance (ρ TOA ) is measured over the Beijing AERONET site. The selected values of ω o and g are used to retrieve the AOD over regions away from the AERONET site and these values are kept spatially constant for a particular image scene. The values ω o and g do not necessarily accurately represent the actual optical properties of the aerosol, since they are influenced by the simplifications involved in the SARA approach. The SARA algorithm is described in detail in [24,26].
In summary, the SEMARA AOD (Equation (1a,b)) at 550 nm was retrieved for

Methods Used for Statistical Analysis
AERONET Sunphotometers do not provide AOD data at 550 nm; therefore, AOD data were interpolated to 550 nm using the Ångström exponent determined from the AOD at 440 and 675 nm (α 440-675 ) and the AOD at 500 nm. To increase the number of statistical samples and to consider the spatial variability imposed by atmospheric motion, AERONET AOD (at least two measurements must be available) was averaged within ± 60 min of the Landsat 8, MODIS, and VIIRS overpass times over Beijing. The SEMARA Landsat , SEMARA MODIS , and SEMARA VIIRS retrieved AOD were averaged over 3 × 3 pixels (at least 2 out of 9 pixels must be available) centered on the AERONET site. It should be noted that different spatial coverage is used with respect to the spatial resolution of the satellite data to investigate the effect of spatial coverage during data sampling on the aerosol retrieval algorithm. For example, spatial coverages used in data sampling for Landsat 8, MODIS, and VIIRS are 0.0081 km 2 , 2.25 km 2 , and 5.06 km 2 , respectively. SEMARA AOD retrievals were validated against AERONET AOD measurements obtained at two sites (Beijing_RADI and Beijing_CAMS) located in Beijing, a city with mixed bright urban land surfaces where the AOD is greatly influenced by severe dust storms and haze episodes. The slope (Equation (8)) and intercept (Equation (9)) between the SEMARA AOD and AERONET AOD were calculated using the reduced major axis (RMA) method, which can simultaneously account for errors in both x and y variables [10,19,37]: where |r| is the absolute value of the Pearson correlation coefficient (r), and S XY and S XX can be calculated using Equations (8a) and (8b), respectively.
where X = mean of X (AERONET) variable, and Y = mean of Y (SEMARA) variable.
The quality and errors of the SEMARA AOD retrievals are reported using the Expected Error (EE, Equation (10)), the root mean squared error (RMSE, Equation (11)), and the relative mean bias (RMB, Equation (12)). EE is the envelope encompassing both absolute and relative error and contains 68% (±1 σ) of the good quality retrievals based on global validation [9,38,39]. Here, the EE for the MODIS Dark Target algorithm at 3 km resolution [8,9,12,[40][41][42] is used to report the good quality SEMARA AOD retrievals.
The upper and lower EE boundaries are calculated using Equations (10a) and (10b), respectively.
RMSE = 0 represents the collocated points on the 1:1 (x = y) line, and RMSE > 0 represents the collocated points scattered away from the 1:1 line.
where RMB > 0.0 and RMB < 0.0 represent over-and under-estimation of AOD retrievals, respectively. The performance of the SEMARA AOD retrievals was also evaluated for different aerosol loadings (low: AOD < 0.2; moderate: 0.2 < AOD < 0.4; high: AOD > 0.4) as well as for different land cover types. Land surface types are classified into three classes using NDVI, i.e., bright surfaces (NDVI ≤ 0.2), sparse vegetation (0.2 < NDVI ≤ 0.4), and moderate vegetation (0.4 < NDVI ≤ 0.6).  Figure 1c) against AERONET AOD measurements at two AERONET sites (Beijing_CAMS and Beijing_RADI) located in Beijing, in an area with mixed and bright urban surfaces. In Figure 1, the solid and dotted black lines represent the 1:1 (y = x) and EE lines, respectively, and the solid red line represents the regression line calculated using the RMA regression (Equations (9) and (10)). In total, 126, 248, and 244 collocated AOD data pairs were available for Landsat-AERONET, MODIS-AERONET, and VIIRS-AERONET, respectively. The results show that SEMARA Landsat , SEMARA MODIS , and SEMARA VIIRS AOD retrievals are well-correlated with AERONET AOD measurements over the urban surfaces of Beijing as indicated by the high values of the Pearson's correlation coefficients, i.e., r = 0.988, 0.980, 0.981, and the slopes of 0.99, 1.01, and 1.0, respectively. For all satellite data, with different spatial resolutions, the collocated points were scattered around the 1:1 line as also indicated by the small values of RMSE = 0.08 for SEMARA Landsat , 0.09 for SEMARA MODIS , and 0.08 for SEMARA VIIRS . A direct relationship was observed between the spatial resolution/coverage of the AOD retrieval and the RMB, i.e., the finer the spatial resolution, the higher the RMB value. For example, a higher RMB of 4.33% was observed for SEMARA Landsat AOD (spatial sampling window: 3 × 3 = 0.0081 km 2 ) and a lower RMB value (−0.54%) was observed for SEMARA VIIRS AOD (spatial sampling window: 3 × 3 = 5.06 km 2 ). This suggests that the resolution of the satellite sensor plays an important role in aerosol inversion. The overestimation of 4.33% in the SEMARA Landsat AOD retrievals may be attributable to a slight underestimation in the SREM estimated surface reflectance during moderate to low aerosol loadings (AOD < 0.25), compared to very high aerosol loading (AOD > 0.50). This overestimation also resulted in a somewhat higher percentage of SEMARA Landsat -retrieved AOD above the EE, i.e., aEE% = 14.29, than for the AOD Remote Sens. 2022, 14, 373 8 of 14 retrieved using SEMARA MODIS (aEE% = 6.86), and SEMARA VIIRS (aEE% = 6.97). Overall, the SEMARA approach performance over the urban surfaces of Beijing was robust as indicated by the large percentage of good quality AOD retrievals within the EE, i.e., wEE% = 85.71, 91.53, and 90.16 for SEMARA Landsat , SEMARA MODIS , and SEMARA VIIRS , respectively. Figure 1d shows the combined validation of the SEMARA Landsat , SEMARA MODIS , and SEMARA VIIRS AOD retrievals. Overall, the retrieved AOD, with different spatial resolution, was well-correlated with AERONET AOD with r = 0.982, slope = 1.01, intercept = 0.003, RMSE = 0.086, RMB = 1.31%, and 89.81% of the retrievals were within the EE (wEE%). SEMARA was tuned to the AERONET AOD at the Beijing site and the results showed the performance of SEMARA in an urban area with heavy traffic (the ring road and other major roads) and, thus, captured well the variety of aerosol concentrations. These results suggest that high-resolution AOD for multi-resolution satellite remote sensing data can be obtained by application of SEMARA over an urban surface with varying and relatively bright surfaces.

Validation of SEMARA AOD Retrievals
in aerosol inversion. The overestimation of 4.33% in the SEMARALandsat AOD retrievals may be attributable to a slight underestimation in the SREM estimated surface reflectance during moderate to low aerosol loadings (AOD < 0.25), compared to very high aerosol loading (AOD > 0.50). This overestimation also resulted in a somewhat higher percentage of SEMARALandsat-retrieved AOD above the EE, i.e., aEE% = 14.29, than for the AOD retrieved using SEMARAMODIS (aEE% = 6.86), and SEMARAVIIRS (aEE% = 6.97). Overall, the SEMARA approach performance over the urban surfaces of Beijing was robust as indicated by the large percentage of good quality AOD retrievals within the EE, i.e., wEE% = 85.71, 91.53, and 90.16 for SEMARALandsat, SEMARAMODIS, and SEMARAVIIRS, respectively. Figure 1d shows the combined validation of the SEMARALandsat, SEMARAMODIS, and SE-MARAVIIRS AOD retrievals. Overall, the retrieved AOD, with different spatial resolution, was well-correlated with AERONET AOD with r = 0.982, slope = 1.01, intercept = 0.003, RMSE = 0.086, RMB = 1.31%, and 89.81% of the retrievals were within the EE (wEE%). SEMARA was tuned to the AERONET AOD at the Beijing site and the results showed the performance of SEMARA in an urban area with heavy traffic (the ring road and other major roads) and, thus, captured well the variety of aerosol concentrations. These results suggest that high-resolution AOD for multi-resolution satellite remote sensing data can be obtained by application of SEMARA over an urban surface with varying and relatively bright surfaces.

SEMARA Performance for Different Aerosol Loadings
Several studies have shown that the accuracy of the retrieved AOD is affected by the aerosol loading in the atmospheric column [10,13,17,18,26,43]. In this section, we examine the performance of SEMARA retrievals for different aerosol loading (low: AOD ≤ 0.2; Remote Sens. 2022, 14, 373 9 of 14 moderate: 0.2 < AOD ≤ 0.4; high: AOD > 0.4). At low aerosol loading, uncertainties related to the surface reflectance estimation have a large influence on the AOD retrieval due to the relatively large contribution of the surface reflectance to the TOA reflectance. However, at high aerosol loading, aerosol contributes relatively more to the TOA reflectance, and the estimation of the microphysical properties of the aerosol particles greatly affects the AOD retrieval [18,44]. Table 2 shows that, for low aerosol loading, 88% of the combined SEMARA AOD retrievals for all sensors (Landsat 8, MODIS, and VIIRS) were within the EE with small RMSE (~0.05) and RMB (~2.50%). This indicates that in the SEMARA approach, the sensitivity of the AOD retrieval to the accurate estimation of the surface reflectance was effectively reduced by using an accurate AOD observation to estimate the aerosol parameters (ω o and g Equation (1a,b)) used for AOD retrieval at other locations than the reference site. With increased aerosol loading and, thus, a relatively larger contribution of aerosol particles to the TOA reflectance, the performance of SEMARA was also good, with a large fraction of the AOD retrievals within the EE (wEE% = 88-94) and low RMSE (0.07-0.13). These results suggest that the SEMARA retrieval approach is capable of retrieving high AOD values without prior information about the aerosol microphysical properties. Table 2. Performance of the SEMARA AOD retrievals during low, moderate, and high aerosol loading conditions. N represents the total number of collocations; aEE%, bEE%, and wEE% represent the percentage of retrievals above, below, and within the expected error; RMSE represents the root mean squared error and RMB represents relative mean bias.

SEMARA Performance over Different Land Cover Types
The brightness of the land surface has a large influence on the AOD retrieval accuracy. The aerosol signal at the TOA over dark surfaces is relatively stronger than over bright surfaces [9,13,14,[45][46][47]. Therefore, in this section, we examine the sensitivity of the SEMARA-retrieved AOD to the surface brightness, using the NDVI as a proxy, i.e., bright surfaces (NDVI < 0.2), sparse vegetation (0.2 < NDVI < 0.4), and moderate vegetation (0.4 < NDVI < 0.6). The greenness of the surface (NDVI) is larger over Beijing_RADI than over Beijing_CAMS, and the NDVI of Beijing_RADI varies seasonally with higher greenness during the spring and early summer. Figure 2 and Table 3 show that, for all three sensors, a larger number of collocated AOD retrievals was obtained for NDVI ≤ 0.2 than for 0.2 < NDVI ≤ 0.4 and for 0.4 < NDVI ≤ 0.6. The SEMARA approach for Landsat 8 performed better over bright surfaces and sparse vegetation than over moderate vegetation ( Table 3). The relatively large error over moderate vegetation might be due to underestimation in SREM estimated surface reflectance, leading to overestimation in the SEMARA AOD retrievals, as indicated by aEE%~25 and RMB~13.35%. However, only 20 collocated Landsat-AERONET observations are available over these moderately vegetated surfaces.    The overall performance of the combined SEMARA-retrieved AOD for all three sensors (Landsat 8, MODIS, and VIIRS) was good over all land surfaces, with 88-92% of the retrievals within the EE, small RMSE (~0.06-0.10), and RMB (−0.54-7.44%) ( Table 3). These results suggest that the SEMARA approach is capable of retrieving AOD over different types of land surfaces. This study shows the robust performance of the SEMARA AOD retrieval approach over a highly polluted city (Beijing, China), which suggests that it can also be implemented in other geographical areas where a reference AOD value is available.

Conclusions
The SEMARA approach, based on the SREM and SARA methods, is introduced for the retrieval of AOD at 30 m resolution from Landsat 8 OLI, 500 m resolution from MODIS, and 750 m resolution from VIIRS. The SEMARA AOD retrievals were validated against the AOD measured at two AERONET sites with mixed and bright urban surfaces located in Beijing. The results show that the SEMARA-retrieved AOD is well-correlated with the AERONET AOD with high values of Pearson's correlation coefficient (r), low values of the RMSE and RMB, and a large percentage of retrievals within the EE (wEE%). A relatively large RMB value was observed for the SEMARA Landsat AOD when the surface was relatively dark (moderate vegetation), as compared to bright surfaces and sparse vegetation. This might be due to underestimation of the SREM estimated surface reflectance and this needs to be thoroughly investigated in future studies. Overall, the SEMARA approach performed well during different aerosol loading as well as over different land surfaces as characterized by the NDVI values. Therefore, the SEMARA approach is highly recommended for aerosol retrieval over mixed surfaces with low to high aerosol loading, using multi-resolution remote sensing data. SEMARA requires an accurate reference AOD value, which could be obtained from ground-based observations or satellite retrievals. The added value to satellite retrieved AOD values is the temporal extension both before and after the satellite overpass.