Window-Based Filtering Aerosol Retrieval Algorithm of Fine-Scale Remote Sensing Images: A Case Using Sentinel-2 Data in Beijing Region

Abstract

The satellite-based Aerosol Optical Depth (AOD) retrieval algorithms are generally needed to construct Land Surface Reflectance (LSR) database. However, errors are unavoidable due to the surface complexity, especially for the short observation period and high-resolution images, such as Sentinel-2 Multi-Spectral Instrument (MSI) data. To address this, reference day images are used instead of the LSR database. The surface is assumed to be Lambertian; however, the fact is that not all pixels meet it well. Therefore, we proposed a window-based AOD retrieval algorithm, which can ignore the unreliable/non-Lambertian pixels in a retrieval window based on two main filtering processes. Finally, using Sentinel-2 Band 1 (60 m), the AODs (120 m) of 134 reference images to 43 reference images were retrieved by this algorithm from 2017 to 2021 in Beijing region, China. The results show that the retrieved AOD with the proposed algorithm exhibits good agreement with the ground-based measured AOD (R > 0.97). The high-resolution AOD presents comparable spatial distributions to the Multiangle Implementation of Atmospheric Correction (MAIAC) algorithm AOD (1 km) products. Moreover, the very little noise and very high spatial continuity of retrieval AOD imply that this algorithm could be ported to other algorithms as part of improving AOD quality.

1. Introduction

Aerosols are colloids of solid particles or liquid droplets suspended in the atmosphere, and typical radii range from 0.001 to 10 μm [1]. They, directly and indirectly, affect radiation forcing, global climate, atmospheric environment and human health [2,3,4].

According to origins, tropospheric aerosols could be divided into two categories: natural sources (such as desert dust) and anthropogenic sources (such as industrial and combustion processes) [5]. Aerosols emitted by the latter usually present high spatial-temporal variability in aerosol optical and radiative properties due to their short lifetime, which also play an important role in many aspects [2]. Therefore, high spatial resolution/fine-scale regional and global aerosol observations are the essential basis of subsequent research.

Aerosol Optical Depth (AOD) is a key parameter that reflects aerosol optical characteristics, representing the extinction of solar radiation caused by aerosol and integrated into the whole atmospheric column and reflecting aerosol content to a certain extent [6]. It could be accessed by one method, i.e., ground-based remote sensing AOD retrieval, with high precise and temporal resolution in a point. However, present observation networks of this method such as AERONET (AErosol RObotic NETwork) and SONET (Sun–Sky Radiometer Observation Network) are too sparse to monitor AOD on a large spatial scale [7]. As an alternative and a deficiency, another method for this, i.e., satellite-based remote sensing AOD retrieval, has been developed and applied to regional and global aerosol studies wildly [8,9,10]. However, due to sensor calibration errors, cloud/cloud shadow contamination, aerosol model classification, and surface reflectance characterization, the satellite-based retrieval AOD has larger uncertainties, especially for fine-scale AOD (less than 1 km) retrieval [11,12].

Until now, researchers have made extensive efforts to retrieve fine-scale and continuous AOD spatial distribution by using the Top-of-Atmosphere (TOA) reflectance of satellite-based observations [4,13]. Most of them first constructed a Land Surface Reflectance (LSR) database, and then it was applied to retrieve AOD [7,14,15,16]. For example, given a satellite TOA image, the LSR could be estimated by two methods according to the density of area vegetation coverage. For densely vegetated areas (area (1)), the visible channel LSR can be estimated by short-wave infrared (SWIR) channels TOA [17]; for high LSR areas with sparse or little vegetation coverage (area (2)), its LSR could be represented by minimum reflectance or the second lowest surface reflectance of a precalculated monthly LSR products [14,18]. Based on these, Wei et al. [14] constructed Landsat LSR database and successfully retrieved 30 m spatial resolution AOD. Validation results showed that the Landsat AOD retrievals were well correlated with the AERONET AOD measurements (R² = 0.93) and better than MODIS MOD04 AOD products. Furthermore, sparsely non-vegetated area in area (2) was independently classified as area (3) by Lin et al. [19] and considered bi-directional reflectance effects. The retrieved AOD with Landsat LSR database constructed by considering these three coverages exhibited good agreements with the ground-based measured AOD (R² = 0.92), and successfully identified polluted sources in two megacities in China. In all these studies, estimating LSR is essential for aerosol retrieval because its 0.01 estimation errors can lead to ~0.1 errors in AOD retrieval [17]. However, for high spatial resolution, on the one hand, few LSR products could reach the meter level, for example, Sentinel-2 Band 2 was 10 m. Despite this satisfaction, the differences in spectral channels and spectral response functions among different sensors may result in inevitable errors in the retrieval algorithm [20]. On the other hand, the presence of retained clouds, topographic shadow, changing surface (vegetation cover areas) and high reflectance surface (urban building areas) make the LSR construction without LSR products difficult and complex [21], especially for Sentinel-2 with a short observation period.

To address the issues above, in this study, we developed a simple, quick and accurate algorithm for retrieving fine-scale AOD without dependency on other LSR products or LSR construction, and it was applied to Sentinel-2 images in Beijing region, China.

2. Study Area and Datasets

2.1. Study Area

As the capital of China, Beijing (Figure 1) is located in the northern part of the North China Plain and is a typically bright urban area where four ground-based AERONET sites (Beijing, Beijing_PKU, Beijing-CAMS and Beijing_RADI) are located. The study area is the footprint of Sentinel-2 image with the highest overlap with Beijing. Therefore, it includes an extra XiangHe site, with a total of five AERONET sites (for more information, refer to

Table 1. Information about AERONET sites used in this study.

No.	Name	Longitude (°E)	Latitude (°N)	City	Level
1	Beijing	116.38	39.98	Beijing	1.5
2	Beijing_RADI	116.38	40.00	Beijing	1.5
3	Beijing_PKU	116.31	39.99	Beijing	1.5
4	Beijing-CAMS	116.32	39.93	Beijing	1.5
5	xianghe	116.96	39.75	Langfang	1.5

).

2.2. Datasets

2.2.1. Sentinel-2 Images

Sentinel-2 is an Earth observation mission under the Copernicus program of the European Space Agency that focuses on observations of the Earth’s surface to provide relevant telemetry services, such as forest monitoring, land cover change detection, and natural hazard management [22]. It consists of two polar-orbiting satellites, Sentinel-2A (launched on 23 June 2015) and Sentinel-2B (launched on 7 March 2017), each carrying a multispectral instrument (MSI) that acquires 13 spectral bands from VIS to SWIR in a 290 km-wide orbit. Each Sentinel-2 satellite has a revisit time of ten days, which is reduced to five days with the combination of the two satellites. The L1C (Level-1C) data available to the public are geometrically refined apparent reflectance data, including 13 bands, the first of which is sensitive to aerosol [23]. L2A (Level-2A) is the surface reflectance image corrected by the official tool Sen2Cor, which also classifies features and generates Scene Classification (SCL) maps during processing, providing an index of features, such as clouds, cloud shadows and water bodies, for each image element.

In this study, Band 1 of Sentinel-2 L1C is used for retrieval, and the SCL of L2A is used for cloud and water mask. The study period is from January 2017 January to December 2021. For more information, refer to Table 2.

2.2.2. Ground Measurements

AEronet RObotic NETwork (AERONET) is an automated ground-based and globally distributed aerosol monitoring network [24]. It was often used to characterize atmospheric aerosols and validate satellite-based AOD retrievals [25,26]. Four AERONET sites in Beijing (Beijing, Beijing_PKU, Beijing-CAMS, Beijing_RADI) and one in the study area (Xianghe) were used. For the latest AERONET Version 3 data, due to the data missing level 2.0 (quality-assured), level 1.5 (cloud-screened) with better data availability was selected to evaluate the retrieved AOD in this study. Considering AERONET provides AOD not at 550 nm but at 440 nm, 500 nm and 675 nm [25], the Ångström exponent algorithm was used to interpolate AOD to 550 nm by using other provided wavelengths above [27].

2.2.3. MODIS MAIAC AOD Product

The Multi-Angle Implementation of Atmospheric Correction (MAIAC) of MODIS (MODerate resolution Imaging Spectroradiometer) is an algorithm for retrieving AOD values over bright and dark surfaces at 1 km resolution using Moderate-resolution Imaging Spectroradiometer (MODIS) products, which utilize time series to separate aerosol and land reflection contributions in TOA reflection [28,29,30]. It considers the effects of bidirectional surface reflectivity by using Ross-Thick Li-Sparse (RTLS) BRF model [31]. Comprehensive assessments showed that MAIAC AOD product has a high correlation coefficient with ground-based observations [32,33,34,35]. Thus, this dataset is used to verify the reliability of our AOD retrieval algorithm in the region without AERONET sites (Table 2).

Table 2. Summary of the satellite images and products used in this study.

No.	Name	Level	Band	Time Span	Spatial Resolution (m)	Center Wavelength (nm)	Notes
1	Sentinel-2	1C	1	2017–2021	60	443	134 images used for retrieval, and 43 images used as reference
2	Sentinel-2	1C	2	2017	10	490	Sensitive to Vegetation Aerosol Scattering (blue) [36], and used in discussion section as an example.
3	SCL	2A	/	2017–2021	60	/	Used for removing cloud, water, shadow, etc.
4	MAIAC	/	/	2017–2021	1000	550	Used for validation.

Table 2. Summary of the satellite images and products used in this study.

No.	Name	Level	Band	Time Span	Spatial Resolution (m)	Center Wavelength (nm)	Notes
1	Sentinel-2	1C	1	2017–2021	60	443	134 images used for retrieval, and 43 images used as reference
2	Sentinel-2	1C	2	2017	10	490	Sensitive to Vegetation Aerosol Scattering (blue) [36], and used in discussion section as an example.
3	SCL	2A	/	2017–2021	60	/	Used for removing cloud, water, shadow, etc.
4	MAIAC	/	/	2017–2021	1000	550	Used for validation.

3. AOD Retrieval Algorithm

3.1. Theoretical Basis

The signal received by the sensor mainly comes from both surface radiation and atmospheric radiation. Assuming a Lambertian surface under a plane-parallel atmosphere, the reflectance at the top of the atmosphere (TOA) received by the sensor can be expressed by the following equation [37,38]:

$${\rho }^{*}({\theta }_s,{\theta }_v,\varphi )={\rho }^a({\theta }_s,{\theta }_v,\varphi )+\frac{\rho T({\theta }_s)T({\theta }_v)}{1-\rho S}$$

(1)

where ${\theta }_s$ is the solar zenith angle, ${\theta }_v$ is the satellite zenith angle, $\varphi$ is the relative azimuth angle between solar and satellite, ${\rho }^a$ is the atmospheric reflectance, $T({\theta }_s)$ and $T({\theta }_v)$ correspond to downward and upward atmospheric transmittance, respectively, $\rho$ is the Lambertian surface reflectance, and $S$ is the spherical albedo of atmosphere. Both ${\rho }^a$ and $T({\theta }_{s/v})$ are the function of $\tau$ (AOD).

The multiple ground-atmosphere interactions are considered small, which can be ignored by setting the $\rho S$ term approach to zero [39,40]. $T({\theta }_s)T({\theta }_v)$ is written as $T({\theta }_s,{\theta }_v)$ and the Equation (1) can be expressed as the following equation:

$${\rho }^{*}({\theta }_s,{\theta }_v,\varphi )={\rho }^a({\theta }_s,{\theta }_v,\varphi )+\rho T({\theta }_s,{\theta }_v)$$

(2)

Once ${\rho }^a$ and $T$ are determined in Equation (2), ${\rho }^{*}$ can be viewed as a linear expression of $\rho$. Furthermore, let us consider two ${\rho }^{*}$ of the same surface but different observation time $t_1$ and $t_2$, and assume that their $\rho$ remains constant in the $\mathsf{\Delta }t$ ($\mathsf{\Delta }t=t_2-t_1$) period. Based on Equation (2), we can obtain the following linear relationship:

$${\rho }^{*}(t_2)=a\times {\rho }^{*}(t_1)+b$$

(3)

where $a=\frac{T(t_2)}{T(t_1)}=f_a({\tau }_1,{\tau }_2){\tau }_i$ and $b={\rho }^a(t_2)-{\rho }^a(t_1)\frac{T(t_2)}{T(t_1)}=f_b({\tau }_1,{\tau }_2)$, and is the AOD of $t_i$ ($i$ = 1, 2). Both $a$ and $b$ are the function of $\tau$.

Inspired by the Contrast Reduction method by Tanre et al. [37], we assume the ${\tau }_1$ (corresponding to $t_1$) could be provided by ground-based measurements rather than being assumed as zero [41], and finally, ${\tau }_2$ is the only unknown variable in Equation (3), which is the basis of retrieving AOD. Moreover, to fit the linear coefficients ($a$ and $b$), at least two pairs of points (one pair of points generate one Equation (3)) are needed. Assuming that the aerosol properties do not vary considerably within a suitable size window, meaning the linear coefficients for all pixels in this window are close, this window can meet at least two pairs of points well. In order to make the linear coefficients in a window vary as little as possible so that they can be fitted stably, we let ${\tau }_1$ and $\mathsf{\Delta }t$ be as small as possible. A small $\tau$ (such as less than 0.1) usually means small variations in $\tau$ as well as $T(t_1)$ and ${\rho }^a(t_1)$ over a region (such as Sentinel-2 footprint), and a small time interval (such as less than 30 days) between $t_1$ and $t_2$ is set to avoid large surface changes. In this condition, in spite of the existence of variations in $T(t_2)$ and ${\rho }^a(t_2)$, and of little surface changes, all of these in a small window (less than 1 km) are so small that linear fitting coefficients can still be calculated steadily, which can be verified by the actual images easily. And linear relationship in a window is the basis of filtering.

3.2. AOD Retrieval

Instead of constructing Land Surface Reflectance (LSR) datasets, the $\rho$ in Equation (2) is obtained directly from a reference image after atmosphere correction. This LSR may bring in more uncertainties than carefully constructed LSR datasets due to the impacts of retained clouds, topographic shadow, changing of surface features and even non-Lambert surface on AOD retrieval (especially for urban regions) [21,42].

Therefore, here, we attempt to propose a window-based retrieval algorithm, which can avoid the above impacts to the maximum extent possible and retrieve fine-scale and continuous AOD. Besides the assumptions in Section 3.1, however, further assumptions are needed for the proposed algorithm: (i) the surface changes little for a short time interval, less than a month for this study; (ii) the aerosol loading can be considered smooth over a small block, less than 1 km for this study. In general, the surface properties change rapidly in space and rather slowly in time [43], and the shorter the time interval, the less the changes in surface properties. Considering the temporal resolution (about 5 days) of Sentinel-2 images and the prevalence of clouds in the Beijing area, the time interval of the first assumption cannot be infinitely small. And for Beijing regions where vegetations are scarce, the first assumption (a month’s time interval) is reasonable for most periods except for the rapid vegetation growth and withering period. A similar compromise can be found in the aerosol retrieval over the land of NOAA (National Oceanic and Atmospheric Administration) AVHRR (Advanced Very High Resolution Radiometer) [44]. Excluding the impacts of the aerosol emission sources or undetected clouds/snow, the AOD (within 1 km) generally meets the second assumption with limited variance over a small block. the smaller the block, the higher the probability to meet this assumption [45]. Under these assumptions, the detailed AOD retrieval processes are described as follows.

Given $t_1$ (hereafter called reference day, ${\tau }_1$ is known) image (hereafter called reference image) and $t_2$ (hereafter called retrieval day, ${\tau }_2$ is unknown) image (hereafter called retrieval image), there are two windows (retrieval windows) in both images corresponding to the same surface extent, such as Figure 2a,b, with a central red $2\times 2$ window. Each retrieval window is contained within a larger auxiliary window, such as a blue $6\times 6$ window in Figure 2. Then, the AOD retrieval can be summarized in three steps: step (1): filtering; step (2): retrieval; step (3): interpolating.

Notably, it may bring great errors if we perform pixel-by-pixel retrieval because the $\rho$ from the reference image contains uncertainties due to the reasons mentioned before. Therefore, we have designed some steps to address it to a certain extent, which is similar to a simple average resampling process, but the average values of a resampling window (corresponding to the retrieval window of this study) are computed by partial pixels. In a resampling window, it makes full use of the surrounding pixels from itself and pixels from the reference image (step (1)) to remove outlier pixels. All outlier pixels are ignored when computing the average value in both retrieval image and reference images, so we accessed two resampled images. They are more in line with the hypothesis (such as Lambertian) as the cost of degraded resolution. The retrievals are then performed in these two resampled images (step (2)). In the end, for very few windows (its pixels are all outliers), a simple interpolate method is used to fill these missing AODs (step (3)). The details are as follows.

Figure 2. Window-based retrieval process. (a) reference image block; (b) retrieval image block; (c–e) three sub-steps; (f–h) three cases. (2 × 2 red window is named the retrieval window, 6 × 6 blue window is named the auxiliary window).

Figure 2. Window-based retrieval process. (a) reference image block; (b) retrieval image block; (c–e) three sub-steps; (f–h) three cases. (2 × 2 red window is named the retrieval window, 6 × 6 blue window is named the auxiliary window).

3.2.1. Filtering

Not all pixels in the retrieval window take part in the AOD retrieval, and outliers that are not dependable for retrieval due to uncertainties as mentioned before need to be removed. Note that the AOD value of the retrieval window is only determined by the pixels within the retrieval window (red pixels). On the other hand, the pixels within the auxiliary window (blue pixels) only serve as references to help eliminate problematic pixels or outliers within the retrieval window. Therefore, this step, referred to as filtering, consists of three sub-steps (as illustrated in Figure 2c–e) that aim to scrutinize and remove outlier pixels.

In the first sub-step, all pixels are masked according to the SCL maps of Sentinel-2 L2A, including cloud, shadow, snow and water.

Furthermore, the retrieval window may exist in three kinds of undependable pixels, including (1) residual clouds, water and shadow; (2) very bright surfaces which suffer from Bidirectional Reflectance Distribution Function (BRDF) effects and are inconsistent with the Lambertian assumption; and (3) changed surfaces which are inconsistent with the surface constant assumption. Based on Section 3.1, the auxiliary window pixels between reference and retrieval day possess a linear relationship under the Lambertian and surface constant assumptions. In contrast, a few pixels deviating significantly from this linear relationship may not meet the above assumptions and need to be filtered out. Therefore, the linear relationship driven by auxiliary window pixels can serve as a filtering sub-step to remove unreliable pixels in a retrieval window. To accomplish this, the variable $f_{diff}$ is first defined as:

$$f_{diff}(i,j)={\rho }^{*}{}_{ij}(t_2)-(a\times {\rho }^{*}{}_{ij}(t_1)+b)$$

(4)

where $(i,j)$ is the index of auxiliary window, $f_{diff}$ is named as difference function, and the definitions of $a$ and $b$ reference Equation (3). Then, an interval based on $f_{diff}$ values is expressed as:

$$A=[\mu -n\sigma ,\mu +n\sigma ]$$

(5)

where $\mu$ is the mean value of $f_{diff}$ values, and their standard deviation is $\sigma$. $n$ controls the size of this interval. Finally, for any pixel $(i,j)$ in a retrieval window, it can be divided into inlier (retained in the subsequent steps) or outlier (removed in the subsequent steps) as follows:

$$\left\{\begin{array}{l}inlier,f_{diff}(i,j)\in A\hfill \\ outlier,f_{diff}(i,j)\notin A\hfill \end{array}\right.$$

(6)

The last but the most important sub-step was designed to exclude pixel pairs with very large or small differences between reference and retrieval days [21]. It can remove more than three undependable kinds of pixels existing on a larger scale (more than a retrieval window but less than an auxiliary window), such as building roofs and small areas of water. In this case, the fitting coefficients itself in the second sub-step may exist as deviations, and the last filtering sub-step serves as an excellent remedy. This sub-step is the same as the second one, but the difference function ($f_{diff}$) is defined as follows:

$$f_{diff}(i,j)={\rho }^{*}{}_{ij}(t_2)-{\rho }^{*}{}_{ij}(t_1)$$

(7)

Here, to differentiate, $n$ for sub-step (2) and (3) is defined as $n_1$ and $n_2$. Generally, the noise in retrieval AOD will be suppressed by the second filtering sub-step using a linear relationship, and the AOD continuity will be enhanced by the third filtering sub-step based on the differences between reference and retrieval images. From these steps, all the masked pixels and outliers will be removed in the next retrieval step. Now, three cases may happen, corresponding to Figure 2f–h: (1) all points are inliers; (2) some are inliers and some are outliers; (3) all points are outliers. Among these three cases, (1) and (2) will continue with step (2), but (3) will execute step (3) after all retrieval windows for executable step (2) are complete.

3.2.2. Retrieval

Subsequently, the mean value of unmasked inliers that were preserved eventually represents the reflectance of the entire retrieval window. Since the ${\tau }_1$ is known, the ${\tau }_2$ is easily calculated based on a Lookup table (LUT) pre-constructed by 6SV (Vector version of Second Simulation of a Satellite Signal in the Solar Spectrum) [38,46] with a series of ${\tau }_2$ values (0.01, 0.20, 0.25, 0.50, 1.00, 1.50, 2.00) using Equation (3). Given a specific ${\tau }_2$ value (such as 0.25), the $a$ and $b$ can be determined in Equation (3). As the ${\rho }^{*}(t_1)$ is known, we can obtain the one-to-one correspondence between ${\tau }_2$ values and ${\rho }^{*}(t_2)$ values. Furthermore, a quadratic function has been used to fit this functional relationship, ${\tau }_2=F({\rho }^{*}(t_2))$. Finally, by bringing the real ${\rho }^{*}(t_2)$ value into this function ($F$), the ${\tau }_2$ value can be easily calculated.

3.2.3. Interpolating

When resampling an image, interpolating methods are wildly used to preserve image quality [47]. Among these methods, the nearest neighbor is the simplest one. For this study, it is enough to fill the retrieval windows that are not available for step (2) because neighbor pixels have a great similarity at high spatial resolutions. Given a pixel that needs to be interpolated, first, we find its nearest surrounding pixels, and then compute their mean value as the interpolate value. If all the nearest surrounding pixels need to be interpolated as well, we expand the search distance (less than 1 km) until we find a valid AOD value.

In summary, the novel algorithm we developed takes full advantage of the surface reflectance provided by the reference day image, but it is selective, and the selection is based on the linear relationships and the differences that exist between the two images. It is worth noting that, in this algorithm, the sizes of retrieval and auxiliary windows, as well as the search pixel radius of the interpolate, are adjustable. In this study, we set them to 120 m, 900 m, and 900 m, respectively. $n_1$ in Equation (5) is set as 1, and $n_2$ in Equation (5) is set as 0.5.

3.3. Aerosol Model

The aerosol model is an important parameter that affects the precision of AOD retrieval and must be considered carefully [15]. Referring to most studies, the particle size distribution in this study is expressed as a mixture of fine and coarse Log-Norm distribution [48,49]. In the actual atmosphere, the optical properties of aerosols change with variations in their concentration and size distribution, so a proven reasonable hypothesis that both size distribution and refractive indices are functions of AOD was employed here [50]. Finally, a dynamic aerosol model was built by statistical analysis of aerosol optical properties in different AOD ranges [51].

Table 3. Optical Properties of the Aerosol Models Used for Lookup Table.

Model	${\mathit{r}}_{\mathit{v}},\mathsf{\mu }\mathbf{m}$	$\mathit{\sigma }$	${\mathit{V}}_{\mathbf{0}},{\mathsf{\mu }\mathbf{m}}^3{/\mathsf{\mu }\mathbf{m}}^2$	Refractive Index: k
fine	$0.064\tau +0.148$	$0.565$	$0.145\tau +0.001$	$1.506+(-0.003\tau +0.01)i$
coarse	$-0.193\tau +3.285$	$0.625$	$0.085\tau +0.065$

displays the Log-Norm size distribution parameters (volume concentration $V_0$, volume median radius $r_v$, and standard deviation $\sigma$) and the refractive indices derived from Beijing_RADI site which has more data records.

4. Results

4.1. Comparison with AERONET

A total of 134 images for retrieval day (AOD > 0.1) and 43 images for reference day (AOD < 0.1) were obtained, and the proposed algorithm was executed on these retrieval images using the closest reference images (<30 days). The average values of retrieval AODs (120 m) in 300 m spatial range were validated with AERONET AOD measurements from five AERONET sites with a 30 min time interval. Figure 3a–e presents the results of the five sites, and Figure 3f is the result of all five sites. Considering the seasonal variations in aerosol characteristics [52], four colors are used to represent spring (March, April, May), summer (June, July, August), autumn (September, October, November), and winter (December, January, February), respectively. One can find the most AOD results in autumn, which is less affected by clouds and snow. Two commonly used statistical indicators are computed here: (1) the correlation coefficient (R) represents the agreement between satellite retrieved AOD and AERONET AOD; (2) the root mean square error (RMSE) is used to measure the differences between satellite retrieved AOD and AERONET AOD. One can note that there are a few summer points in Figure 3, which are caused by the prevalence of clouds and the difficulty in meeting the 30 day interval.

In Figure 3, all site validation results showed very high agreements (R > 0.95, RMSE < 0.1), and their slopes of trendline are close to $y=x$, especially at the Beijing site. These illustrate both the stability of the algorithm (with high R values for all sites) and the rationality of the aerosol model parameters (the fitted line is close to $y=x$). However, the algorithm tends to overestimate the AOD in spring slightly (as most of the red scattered dots representing spring are above the $y=x$ line), suggesting that a 30 day interval may not meet the assumption of constant surface reflectance during the vegetable growth period, or that the total aerosol model cannot represent this season well.

4.2. Evaluations in Spatial Distribution

To validate the performance of the new algorithm in areas without AERONET sites, the wildly evaluated MODIS MAIAC AOD product was used. Here, we evaluated the spatial distribution of high-value AOD (max value more than 1) and low-value AOD (max value less than 0.5) separately, where the latter is more prone to errors due to the dominance of surface signals.

For high-value AOD (Figure 4), overall, the AOD spatial trends of the two algorithms are consistent. the local differences may be attributed to the overpass time difference between MODIS and Sentinel-2 satellites (~5 h), and the gaps represent masked contaminated pixels (clouds, shadows, or waters). Benefiting from the process of filtering, more details of AOD could be observed with less noise in the new algorithm, which is useful for detecting fine aerosol emission sources.

As with high-value AOD, the low-value AOD (Figure 5) keeps consistent with AERONET in spatial well. This illustrates that the images of reference day we selected are reasonable, and also proves that the proposed algorithm could capture both large and subtle differences between images of retrieval days and reference days. It is clear that there is an obvious discontinuity in our inverse AOD result on 1 November 2020, which can be attributed to the image quality but not the retrieval algorithm. To our knowledge, this phenomenon occurs frequently for retrieval days with rapid variations in a short time.

5. Discussion

Under the assumptions of Lambertian surface, this study proposed a window-based (instead of pixel-based) aerosol optical depth (AOD) retrieval algorithm. Benefiting from the filtering process of this algorithm, it can retrieve image AODs based on a reference image. In fact, from one perspective, similar works (using a reference image to retrieve AODs) can date back to the 1980s. In 1988, Tanre et al. [37] proposed the structure function algorithm, which can derive image AODs over the Saharan area by using structure information (based on windows) of a reference image. Later, many researchers put the effort into improving the complexity of structured functions and applying it to urban areas [40,53,54]. However, these works still use all pixels in a window, some of which cannot meet the Lambertian assumptions and may introduce errors in the calculation of the structure function. From another perspective (removing window outliers), under the assumptions of Lambertian surface, Sun et al. [21] acquired Land Surface Reflectance using the clear sky composite technique. For each window, the 30% brightest and 30% darkest pixels were excluded, and the average value of the remaining 40% pixels was used to represent the reflectance of the entire window. This closely resembles the last filtering step proposed in the study. In addition, when evaluating AOD products, a series of filtering steps are carried out to discard AOD outliers within a window caused by factors such as BRDF (Bidirectional Reflectance Distribution Function), mixed pixels and cloud [55], and to improve verification accuracy [56]. Unfortunately, these steps occur in the validation phase after the retrieval.

To address these limitations above, a window-based AOD algorithm was proposed in this study to filter out these pixels. The proposed algorithm achieved high accuracy (R > 0.95 for all AERONET sites) and spatial continuity in the AOD retrieval.

However, this study still suffers some uncertainties that need further improvement. First, one aspect in the filtering process of the proposed algorithm that should be figured out is the pixel utilization in a retrieval window. To explore this, we implemented an algorithm in Sentinel-2 Band 2 (10 m spatial resolution). Its retrieval window was set at $10\times 10$, so the AOD has a 100 m spatial resolution with the same spatial distribution as 120 m (Figure 4). Due to a $10\times 10$ retrieval window, the pixel number (100 for Band 2) is much larger than that in a $2\times 2$ retrieval window (4 for Band 1), which allows the number distribution map of inliers (Figure 6) in the retrieval window to have more levels. In Figure 6, we selected two regions, where the red region represents vegetable region, and the blue one region is urban region. In region one, the number of buildings, waters and mountain shadows is low, but the number of vegetation areas is high. In region two, the number is overall lower than that in region one due to less vegetation. In comparison, the very bright buildings have lower numbers near the waters. Moreover, in both region one and region two, the retrieval windows without inlier after filtering (the number is 0 and the white color is not interpolated in the AOD map), where all pixel brightness is very high or low compared to the pixels in auxiliary window, occupy a very small percentage, so the interpolated and non-interpolated images are very similar. Based on the above analysis, a potential problem is that the outliers cannot be removed if all pixels in the auxiliary window have very high or low brightness, which may lead to spatial discontinuity on a larger region scale (for example, more than 1 km in this study). In this case, a larger auxiliary window is needed; however, a steady linear relationship decreases as the window increases. Therefore, a suitable auxiliary window size should be identified first before applying this algorithm to other regions.

The second issue is the rationality of the assumption that the surface changes little in one month. In order to validate it, we have selected one retrieval image (24 October 2020) to study the influences as the time interval (between the retrieval image and the reference image) increases. In the study period, this retrieval image has a lot of reference images, and the time intervals are 3 days, 5 days, 13 days, 30 days, 40 days, and 66 days, respectively. Figure 7 shows the retrieval results based on these reference images. Note that the incomplete results for some days (3 days, 13 days, and 66 days) are attributed to the absence of data in the reference or retrieval images due to the orbital design of the satellite. Here, the retrieval result based on the closest reference image (3 days) is seen as the most accurate and is used for comparison with those based on other reference images, as shown in Figure 8. Figure 7 and Figure 8 suggest that the retrieval results are very close until the time interval reaches 30 days. The AOD distribution maps of 30 days and 40 days still resemble that of 3 days, but their scattering points are more discrete. This spatial similarity does not last when the time interval reaches 66 days, whose AOD distribution map is significantly different from 3 days. In general, for this study and this retrieval day, if the time interval is less than 40 days, the influences of surface changes brought to the final AOD results are acceptable.

Another issue is the discontinuity brought about by the remote sensing image itself, particularly in the case of Sentinel-2 Band 1. Although the proposed algorithm can effectively filter out outliers caused by surface features, it cannot address this particular problem. As such, a preprocessing step to mitigate this issue is necessary and is recommended for future work.

Finally, the AOD results during the summer season, when vegetation undergoes rapid changes, were impacted by cloud cover, resulting in fewer retrievals. It remains to be seen whether the proposed algorithm can effectively filter out pixels with significant differences between retrieval and reference images caused by vegetation growth. Further validation is required in future work.

6. Conclusions

In this study, we proposed an aerosol retrieval algorithm using two-day images including a reference day and a retrieval day instead of constructing an LSR database. Considering subtle LSR changes and the BRDF effects, we designed a filtering process in this algorithm to remove unreliable pixels in a retrieval window, including (1) pixels marked cloud, water, shadow, and snow label in SCL; (2) pixels with high or low difference between the two-day images above; and (3) pixels that deviate significantly from the linear relationship between the two-day images. The retrieval AODs (120 m) in Beijing using Sentinel-2 Band 1 exhibit high agreements with AERONET sites (R > 0.95) and common spatial distributions with MODIS MAIAC AOD. Benefiting from the filtering process and high-resolution images, the retrieval AODs show continuous and detailed spatial distributions relative to MAIAC AODs during low AOD and high AOD days.

Some uncertainties still exist and need to be further improved in the proposed algorithm. Firstly, in the filtering process, an auxiliary window is needed, and it is necessary to determine a suitable window size for different regions according to the land surface. Secondly, the errors from the image itself need to be reduced before the retrieval algorithm. Lastly, the AODs in summer in this study are too few to validate the filtering effects of the algorithm for vegetation growth pixels.

Finally, it should be noted that this algorithm is designed for retrieving AODs using a reference image, although it could also be used in algorithms based on an LSR database to obtain greater precision.