An Improved Method for Retrieving Aerosol Optical Depth Using Gaofen-1 WFV Camera Data

: The four wide-ﬁeld-of-view (WFV) cameras aboard the GaoFen-1 (GF-1) satellite launched by China in April 2013 have been applied to the studies of the atmospheric environment. To highlight the advantages of GF-1 data in the atmospheric environment monitoring, an improved deep blue (DB) algorithm using only four bands (visible–near infrared) of GF-1/WFV was adopted to retrieve the aerosol optical depth (AOD) at ~500 m resolution in this paper. An optimal reﬂectivity technique (ORT) method was proposed to construct monthly land surface reﬂectance (LSR) dataset through converting from MODIS LSR product according to the WFV and MODIS spectral response functions to make the relationship more suitable for GF-1/WFV. There is a good spatial coincidence between our retrieved GF-1/WFV AOD results and MODIS/Terra or Himawari-8/AHI AOD products at 550 nm, but GF-1/WFV AOD with higher resolution can better characterized the details of regional pollution. Additionally, our retrieved GF-1/WFV AOD (2016–2019) results showed a good agreement with AERONET ground-based AOD measurements, especially, at low levels of AOD. Based on the same LSR dataset transmitted from 2016–2018 MODIS LSR products, R ORT of 2016–2018 and 2019 GF-1/WFV AOD retrievals can reach up to 0.88 and 0.94, respectively, while both of RMSE ORT are smaller than 0.13. It is indicated that using the ORT method to deal with LSR information can make GF-1/WFV AOD retrieval algorithm more suitable and ﬂexible.


Introduction
A significant portion of aerosols in the atmosphere sources from anthropogenic activities, including industrial processes, fossil fuel combustion, agricultural operation, construction and mining. In recent years, China suffered from frequently severe pollution events, especially in winter, which had strong impact on the air environment, climate change, and public health [1][2][3][4]. However, due to the heterogeneous distribution of sources, short lifetime, and episodic features of emission events, aerosols exhibit high spatiotemporal variability which can hardly be characterized by the sparsely ground-based measurements. Therefore, aerosol optical depth (AOD) retrieved from satellite data has been increasingly used to estimate the surface-level particle concentrations or the air pollution level [5].
Although there are some uncertainties (e.g., differences in the sensor, calibration/ characterization, retrieval algorithm, pixel selection, cloud and other masking) in satellite AOD retrievals, many relevant studies still been continually developed [6][7][8][9][10]. Dark target (DT) and deep blue (DB) algorithms have been successfully applied to the Moderate Resolution Imaging Spectroradiometer (MODIS), Visible Infrared Imaging Radiometer Suite (VIIRS), and Medium Resolution Spectral Imager (MERSI) onboard the Chinese FengYun (FY) satellite series (i.e., FY-3A, FY-3B, FY-3C and FY-3D) and other sensors with medium-low resolution [5,[11][12][13]. MODIS AOD products (MOD/MYD04) provide daily global AOD distributions at 10 and 3 km spatial resolutions since 2000, and they are  Table 1 shows the specific configuration parameters of GF-1/WFV. The data of each GF-1/WFV image includes four files: (1) The GeoTIFF file is used to store image data in a 32-bit integer format; (2) JPG file is used to display thumbnails of images; (3) XML files include data description files, which mainly provide angle (Solar Zenith Angles, Satellite Zenith Angles, Solar Azimuth Angles and Satellite Azimuth Angles) information of image center point, transit time, etc.; (4) The RPB file mainly provides the geometric positioning information of the image, which is used for geometric correction of image.

Other Satellite Products
The DT (AOD_MODT) and DB (AOD_MODB) AOD datasets from Terra/MODIS aerosol products (MOD04), Himawari-8/AHI aerosol product (AOD_AHI) [35], and groundbased AOD observations from AERONET (AOD_AERO) [36,37] were used for the validation of retrieved GF-1/WFV AOD (AOD_WORT) results at 550 nm. Here, the Ångström exponent (AE) relationship was applied for the calculation of AOD_AHI and AOD_AERO at 550 nm. Terra/MODIS surface reflectance product (MOD09A1) [12,14,34], denoted as MODIS LSR. Detailed product information is listed in Table 2.   Table 1 shows the specific configuration parameters of GF-1/WFV. The data of each GF-1/WFV image includes four files: (1) The GeoTIFF file is used to store image data in a 32-bit integer format; (2) JPG file is used to display thumbnails of images; (3) XML files include data description files, which mainly provide angle (Solar Zenith Angles, Satellite Zenith Angles, Solar Azimuth Angles and Satellite Azimuth Angles) information of image center point, transit time, etc.; (4) The RPB file mainly provides the geometric positioning information of the image, which is used for geometric correction of image.

Other Satellite Products
The DT (AOD_MODT) and DB (AOD_MODB) AOD datasets from Terra/MODIS aerosol products (MOD04), Himawari-8/AHI aerosol product (AOD_AHI) [35], and ground-based AOD observations from AERONET (AOD_AERO) [36,37] were used for the validation of retrieved GF-1/WFV AOD (AOD_WORT) results at 550 nm. Here, the Ångström exponent (AE) relationship was applied for the calculation of AOD_AHI and AOD_AERO at 550 nm. Terra/MODIS surface reflectance product (MOD09A1) [12,14,34], denoted as MODIS LSR. Detailed product information is listed in Table 2. Retrieved GF-1/WFV AOD data was matched with the AERONET observations from 2016 to 2019 over BTH region. Within a sampling window of 25 km × 25 km centered in an AERONET site, more than 25% of the total valid pixels are required to calculate a mean value of GF-1/WFV AOD for the accuracy evaluation. And no less than 4 AERONET AODs at 550 nm with a time interval of 1 h centered at Terra overpass time were averaged as the corresponding ground-based measurements for our validation.

Deep Blue (DB) Algorithm
In the case of an infinite uniform Lambertian and homogeneous target of surface reflectance ρ s λ , the reflectance of the top of the atmosphere ρ T λ at a particular wavelength λ can be written as follows [31,38,39]: where ρ a λ is the atmospheric path reflectance, µ s is the cosine of the solar zenith angle (θ s ), µ v is the cosine of the view zenith angle (θ v ), φ is the relative azimuth angle (the difference between the solar azimuth angle and the satellite azimuth angle), T ↑ λ (µ v ) is the total upward transmission in the direction of the satellite field of view, T ↓ λ (µ s ) is the total downward atmospheric transmission, and S λ is the spherical albedo of the atmosphere for illumination from below (atmospheric backscattering ratio). On the left-hand side of Equation (1), ρ T λ can be obtained from satellite observation data. Except for the surface reflectance ρ s λ , each term on the right-hand side of Equation (1) is a function of the aerosol type and loading optical thickness τ.
Using the radiative transport model Second Simulation of Satellite Signal in the Solar Spectrum (6S) [39,40], the relationship between aerosol optical depth and parameters (i.e., S λ ρ atm λ , T ↓ λ (µ s )T ↑ λ (µ v )) was calculated under different atmospheric aerosol modes and observation conditions. A set of sun-satellite geometries, atmospheric parameters, and aerosol information listed in Table 3 was applied to build a look-up table (LUT), which would be used in the AOD retrieval to speed up the calculation processes. The Rayleigh scattering has significant impacts on the visible channels, especially for the blue band with short wavelength [15]. In this study, the method proposed by Levy et al. was used to correct the altitude of atmospheric molecule scattering [5]: where τ λ,z is the Rayleigh optical depth (ROD), λ is the wavelength (µm), Z is the elevation of the surface target obtained from GDEM, and 8.5 km is the exponential scale height of the atmosphere. An overview of the retrieval algorithm is shown in the dataflow diagram in Figure 2. The original GF-1/WFV data was sampled with 32 × 32 pixels to obtain a new image with a resolution of 512 m (~500 m, regarded as 500 m). Then, angle calculation (including view zenith angle, solar zenith angle and relative azimuth angle) and pixel selection were pre-processed before the AOD retrieval. GF-1/WFV surface reflectance converted from the MODIS LSR data and seven standard aerosol models were used to identify the land surface information and aerosol type in the GF-1/WFV AOD calculation, respectively. Here, the above seven aerosol models involved two standard aerosol models (urban and continental models) in 6S [39,40] and five models (dust, smoke-high absorption, smoke-low absorption, urban-clean and urban-polluted models) in the VIIRS over-land aerosol retrieval [41]. Once the calculated spectral reflectance makes the best match with those that are measured, the aerosol model is determined and the AOD value will be reported.

Angle Calculation
For the original GF-1/WFV data, only the angle information in the center of the satellite image is provided. The difference in the observation angle between the image edge pixel and the central pixel can cause a large error to the quantitative application. The variations of apparent reflectance changed with the observation geometry and AOD were simulated by the 6S model, and the results are shown in Figure 3. Figure 3c,d indicate that the change of relative azimuth angle has little influence on the apparent reflectance for a fixed AOD, so that the entire image can utilize the uniform value of relative azimuth angle provided by the XML file. However, based on the simulations shown in Figure 3a,b, the satellite zenith angles of the satellite image need to be calculated by interpolation of camera field angles due to the strong effect of satellite zenith angle on the apparent reflectance. According to Jacobson (2005) [42], the solar zenith angle can be calculated by using the

Angle Calculation
For the original GF-1/WFV data, only the angle information in the center of the satellite image is provided. The difference in the observation angle between the image edge pixel and the central pixel can cause a large error to the quantitative application. The variations of apparent reflectance changed with the observation geometry and AOD were simulated by the 6S model, and the results are shown in Figure 3. Figure 3c,d indicate that the change of relative azimuth angle has little influence on the apparent reflectance for a fixed AOD, so that the entire image can utilize the uniform value of relative azimuth angle provided by the XML file. However, based on the simulations shown in Figure 3a,b, the satellite zenith angles of the satellite image need to be calculated by interpolation of camera field angles due to the strong effect of satellite zenith angle on the apparent reflectance. According to Jacobson (2005) [42], the solar zenith angle can be calculated by using the satellite transit time, longitude and latitude of the corresponding satellite image pixel, as follows: cos(θ s ) = sin(lat) × sin δ + cos(lat) × cos δ × cos t where δ is the declination of the sun, t is the solar hour angle, and lat is the latitude of the pixel. The longitude and latitude of the image elements can be bilinearly interpolated by the longitude and latitude values of four image corners.

Appropriate Pixel Selection
In this paper, some obstructions such as cloud, water, and snow/ice pixels need to be identified before AOD retrieval. The premise of AOD retrieval is the accurate isolation of cloud and clear regions. Due to the limitation of GF-1/WFV band, it is difficult to effectively remove cloud pixels but maintain heavy aerosol ones. Figure 4 shows the GF-1/WFV cloud recognition flow, which is described in detail below.
The multi-band threshold approach of cloud removing, which is combined with visible bands and infrared bands and commonly applied to the MODIS or Landsat data, is not suitable for the GF-1/WFV with only four bands. Bright clouds characterized by high reflectivity in green and red channels, can lead to the signal oversaturation and even an upward trend of reflectance spectrum. Therefore, a combination approach of multi-band threshold method and reflectance variation characteristics with the increasing wavelength can be used for the cloud identification of GF-1/WFV. Due to differences in the microphysical and chemical properties between cloud and aerosol, cloud exhibit a more complex spatial behavior than aerosol [43]. In our study, interesting regions, including thick cloud, thin cirrus cloud and clear sky, were extracted

Appropriate Pixel Selection
In this paper, some obstructions such as cloud, water, and snow/ice pixels need to be identified before AOD retrieval. The premise of AOD retrieval is the accurate isolation of cloud and clear regions. Due to the limitation of GF-1/WFV band, it is difficult to effectively remove cloud pixels but maintain heavy aerosol ones. Figure 4 shows the GF-1/WFV cloud recognition flow, which is described in detail below.
The multi-band threshold approach of cloud removing, which is combined with visible bands and infrared bands and commonly applied to the MODIS or Landsat data, is not suitable for the GF-1/WFV with only four bands. Bright clouds characterized by high reflectivity in green and red channels, can lead to the signal oversaturation and even an upward trend of reflectance spectrum. Therefore, a combination approach of multi-band threshold method and reflectance variation characteristics with the increasing wavelength can be used for the cloud identification of GF-1/WFV. Due to differences in the microphysical and chemical properties between cloud and aerosol, cloud exhibit a more complex spatial behavior than aerosol [43]. In our study, interesting regions, including thick cloud, thin cirrus cloud and clear sky, were extracted by a supervised classification method [44]. Due to the difference in the apparent reflectance of the image under different backgrounds, histogram statistics were conducted for exploring the pixel distribution to determine the threshold for dividing clouds and non-clouds. The spatial variation of the apparent reflectance at visible wavelength is relatively uniform for clear sky pixels, while that of cloud pixels is extremely mutational. And with the comparison of clear-sky condition, the amplitude of spatial variation of apparent reflectance can be further deduced for heavy aerosol-loaded pixels due to the stronger scattering and absorption led by aerosol particles. Therefore, the standard deviation of band 1 (0.484 µm) TOA reflectance within a window of 3 × 3 pixels (3 × 3-STD) is higher for cloud pixels than that for non-cloud pixels [34]. To eliminate the misclassification of cloud pixels, isolated cloud pixels were regrouped into cloud-free type if the proportion of cloud pixel number within a window of 3 × 3 pixels is less than 25%. A buffer zone with a distance of 3-pixel size was used for expanding the initially recognized cloud regions to guarantee that the cloud pixels in the satellite image can be removed completely. flectance can be further deduced for heavy aerosol-loaded pixels due to the stronger scattering and absorption led by aerosol particles. Therefore, the standard deviation of band 1 (0.484 μm) TOA reflectance within a window of 3 × 3 pixels (3 × 3-STD) is higher for cloud pixels than that for non-cloud pixels [34]. To eliminate the misclassification of cloud pixels, isolated cloud pixels were regrouped into cloud-free type if the proportion of cloud pixel number within a window of 3 × 3 pixels is less than 25%. A buffer zone with a distance of 3-pixel size was used for expanding the initially recognized cloud regions to guarantee that the cloud pixels in the satellite image can be removed completely. Because the algorithm is sensitive to snow/ice and water, the normalized difference vegetation index (NDVI) is a simple graphical indicator that can be further used to remove those non-inversion pixels. Cloud, water, and snow/ice pixels can be also identified based on a rule of NDVI < 0.

Reflectance Relationship Transformation
In this paper, a new aerosol retrieval algorithm was proposed to retrieve GF-1/WFV AOD based on the MODIS LSR data. As a key process, the transformational relationship between the surface reflectance of GF-1/WFV and MODIS was implemented through two main steps: the building of monthly MODIS LSR dataset and the establishment of LSR transfer coefficients.
MODIS LSR data were obtained from eight-day synthesized MODIS 500 m surface reflectance product after atmospheric corrections for gaseous scattering and absorption, aerosol scattering and absorption, cirrus contamination, BRDF coupling, and the adja- Because the algorithm is sensitive to snow/ice and water, the normalized difference vegetation index (NDVI) is a simple graphical indicator that can be further used to remove those non-inversion pixels. Cloud, water, and snow/ice pixels can be also identified based on a rule of NDVI < 0.

Reflectance Relationship Transformation
In this paper, a new aerosol retrieval algorithm was proposed to retrieve GF-1/WFV AOD based on the MODIS LSR data. As a key process, the transformational relationship between the surface reflectance of GF-1/WFV and MODIS was implemented through two Remote Sens. 2021, 13, 280 9 of 19 main steps: the building of monthly MODIS LSR dataset and the establishment of LSR transfer coefficients.
MODIS LSR data were obtained from eight-day synthesized MODIS 500 m surface reflectance product after atmospheric corrections for gaseous scattering and absorption, aerosol scattering and absorption, cirrus contamination, BRDF coupling, and the adjacency effect [40]. In the first step, the dataset of 500 m monthly MODIS LSR in blue wavelength was constructed by using the traditional minimum reflectivity technique (MRT) and the optimal reflectivity technique (ORT) method. The MRT method considers a certain number of images at different times, and seeks the minimum value in the same position pixel as the target value [25,45]. Although the monthly LSR data constructed by the MRT method is the LSR minimum of each pixel in one month, the actual land surface reflectance may be underestimated. By comparison, although the ORT method cannot guarantee the existence of every pixel, the generated LSR data is more rigorous to ensure the accuracy of AOD inversion. The determination rules of ORT method are as follows: MOD09A1 data from 2016 to 2018 were selected, and pixels with view zenith angle greater than 30 degree were removed. Each pixel is sorted from low to high values, then set the rules: where (i,j) is the position of a single pixel, L(i,j) is the unsorted LSR value, and I(i,j) is the sorted LSR value. N is the number of primitive values per month for the same location, n is the number of values processed (view zenith angle ≤ 30 • ) per month at the same location. S(i,j) is the final MODIS LSR which can be determined as follows: where M k is an array containing three elements L 3k (i, j), L 3k−1 (i, j) and L 3k−2 (i, j). According to the significant difference between MODIS and GF-1/WFV bands in the wavelength range and spectral response function (SRF) shown in Figure 5, the surface reflectance of GF-1/WFV can be inherited from MODIS LSR, but cannot be replaced directly. Theoretically, the Earth observation of the four WFV sensors aboard on GF-1 satellite is simultaneous, meaning that the observation geometries of four WFVs are not the same. However, our simulated results in Figure 5b indicate that for the blue band, the impact of four different WFV cameras on the apparent reflectance can be ignored for a single observation of GF-1/WFV. Therefore, the spectral differences among the four cameras were not considered into the transformation method of surface reflectance, and a single linear relationship between MODIS and GF-1/WFV surface reflectance was assumed in the second step.
where A and B are two coefficients of the linear formula, and ρ s_MODIS λ and ρ s_WFV λ are surface reflectance values observed by MODIS and WFV, respectively, which can be derived by: where Γ(λ) is the SRF. λ 1 and λ 2 are the upper and lower limits of integration, respectively, which define the band width, and S(λ) is the spectral curve corresponding to central wavelength λ.
After differentiation, Equation (6) can be rewritten as: where ∆λ = 0.0001. In this study, the spectral library established in this paper is obtained by using the remote sensing software ENVI. The spectral data of 50 typical features (construction: 15 types, soils: 20 types, and vegetation: 15 types) in the ENVI spectrum database (ENVI standard spectral library (SLI) files, THOR Metadata Rich Spectral Library (MRSL) output, and Analytical Spectral Devices (ASD) spectrometers output) [34] was selected to calculate the reflectance of ground objects in blue band for both GF-1/WFV and MODIS sensors. As shown in Figure 6 and Table 4, the surface reflectance of GF-1/WFV is higher than that of MODIS, and there is little difference in the transformation parameters of the four WFV cameras. Then, by fitting all the scattered points in Figure 6, the GF-1/WFV blue surface reflectance relationship can be available from Equation (9). In this study, the spectral library established in this paper is obtained by using the remote sensing software ENVI. The spectral data of 50 typical features (construction: 15 types, soils: 20 types, and vegetation: 15 types) in the ENVI spectrum database (ENVI standard spectral library (SLI) files, THOR Metadata Rich Spectral Library (MRSL) output, and Analytical Spectral Devices (ASD) spectrometers output) [34] was selected to calculate the reflectance of ground objects in blue band for both GF-1/WFV and MODIS sensors. As shown in Figure 6 and Table 4, the surface reflectance of GF-1/WFV is higher than that of MODIS, and there is little difference in the transformation parameters of the four WFV cameras. Then, by fitting all the scattered points in Figure 6, the GF-1/WFV blue surface reflectance relationship can be available from Equation (9). late the reflectance of ground objects in blue band for both GF-1/WFV and MODIS sensors. As shown in Figure 6 and Table 4, the surface reflectance of GF-1/WFV is higher than that of MODIS, and there is little difference in the transformation parameters of the four WFV cameras. Then, by fitting all the scattered points in Figure 6, the GF-1/WFV blue surface reflectance relationship can be available from Equation (9).

Results and Analysis
Here, ORT LSR results were firstly compared with the MRT LSR data. Then, GF-1/WFV AODs over BTH from 2016 to 2019 were retrieved based on the new algorithm mentioned in Section 3. To assess the accuracy of our retrieved AOD_WORT, 550 nm AODs obtained from AERONET ground-based measurements and different satellite observations were used for the validation.

Land Surface Reflectance
In order to verify the role of LSR data in the GF-1/WFV AOD retrieval, a sensitive analysis was carried out by simulating the apparent reflectance in GF-1/WFV blue band under different AOD and LSR conditions. As shown in Figure 7, the apparent reflectance of GF-1/WFV in blue band is very sensitive to the LSR when the AOD is less than 3. However, with the increase of AOD, the apparent reflectance grows slowly for all LSR conditions. Therefore, an inaccurate LSR data can introduce a large uncertainty into the GF-1/WFV AOD retrieval, and the proper estimation of LSR which can better capture the characteristics of actual surface condition is very important for the real-time observation of GF-1/WFV AOD.

Results and Analysis
Here, ORT LSR results were firstly compared with the MRT LSR data. Then, GF-1/WFV AODs over BTH from 2016 to 2019 were retrieved based on the new algorithm mentioned in Section 3. To assess the accuracy of our retrieved AOD_WORT, 550 nm AODs obtained from AERONET ground-based measurements and different satellite observations were used for the validation.

Land Surface Reflectance
In order to verify the role of LSR data in the GF-1/WFV AOD retrieval, a sensitive analysis was carried out by simulating the apparent reflectance in GF-1/WFV blue band under different AOD and LSR conditions. As shown in Figure 7, the apparent reflectance of GF-1/WFV in blue band is very sensitive to the LSR when the AOD is less than 3. However, with the increase of AOD, the apparent reflectance grows slowly for all LSR conditions. Therefore, an inaccurate LSR data can introduce a large uncertainty into the GF-1/WFV AOD retrieval, and the proper estimation of LSR which can better capture the characteristics of actual surface condition is very important for the real-time observation of GF-1/WFV AOD.
According to the ORT method, it is not difficult to understand that ORT LSR is always higher than MRT LSR. By comparing the monthly ORT LSR results over BHT region with the corresponding MRT ones, as shown in Figure 8, the difference of monthly LSR between MRT and ORT methods is in a range of −0.031 ± 0.014.
BTH vegetation with obvious zonation was predominated with deciduous broadleaved forest, and the vegetation cover changes significantly with the seasons. Due to the growth change of land surface vegetation, the land surface changes obviously in spring and autumn. Additionally, the land surface covered by snow and ice differs significantly from that covered by bare soil or wheat in winter. Therefore, the minimum reflectance may not be a reasonable representation of the actual surface reflectance. Although there are random errors in AOD retrieval by using ORT LSR, these errors can be reduced by taking the average. In summer, the difference between MRT and ORT LSRs is minimal due to the small change of underlying surface (Figure 8).  According to the ORT method, it is not difficult to understand that ORT LSR is always higher than MRT LSR. By comparing the monthly ORT LSR results over BHT region with the corresponding MRT ones, as shown in Figure 8, the difference of monthly LSR between MRT and ORT methods is in a range of −0.031 ± 0.014.

Spatial Distributions of AODs
To illustrate the reliability of our retrieved GF-1/WFV AOD_WORT in terms of spatial distribution, an inter-comparison was performed based on GF-1/WFV, Terra/MODIS, Himawari-8/AHI AODs derived from different retrieval algorithms. Figures 9 and 10 depict the spatial distributions of AOD_WMRT (500 m), AOD_WORT (500 m), AOD_MODT (10 km), AOD_MODB (10 km), and AOD_AHI (0.05 deg) on 12 July 2017 and 1 January 2018, respectively. The heavily polluted areas were mainly concentrated on the southeast part of BTH region, which were obviously different from the northwest mountainous areas. The spatial distribution of our AOD_WORT retrievals is generally consistent with those of other four satellite AODs. The high values of AOD (>0.8) are mainly distributed in the southeast of BTH region, and zones with the low AODs are located in the northwest parts. In the northwestern BTH region with lower AOD values area, AOD_WORT, and AOD_MODB are much smoother than other three AOD results. It is notable that the invalid value ratio of GF-1/WFV AODs is much lower than those of Terra/MODIS, Himawari-8/AHI ones, especially for AOD_MODT and AOD_AHI, which may be caused by the difference in spatial resolution, overpassing time, and retrieval algorithm. Such a situation of a high invalid value ratio was particularly serious for the AOD_MODT over the southeastern BHT region with air pollution, which made the status over this region cannot be captured for the analysis. And some pixels covered with snow and ice in Figures  9b,c and 10b,c were also set as invalid values the northwest mountains. However, for the retrieved GF-1/WFV AOD with a spatial resolution of 500 m, both AOD_WORT and AOD_WMRT can provides more abundant details than MODIS and AHI in the spatial distribution, and can more accurately reflect the AOD spatial variation over the polluted region. BTH vegetation with obvious zonation was predominated with deciduous broadleaved forest, and the vegetation cover changes significantly with the seasons. Due to the growth change of land surface vegetation, the land surface changes obviously in spring and autumn. Additionally, the land surface covered by snow and ice differs significantly from that covered by bare soil or wheat in winter. Therefore, the minimum reflectance may not be a reasonable representation of the actual surface reflectance. Although there are random errors in AOD retrieval by using ORT LSR, these errors can be reduced by taking the average. In summer, the difference between MRT and ORT LSRs is minimal due to the small change of underlying surface (Figure 8).

Spatial Distributions of AODs
To illustrate the reliability of our retrieved GF-1/WFV AOD_WORT in terms of spatial distribution, an inter-comparison was performed based on GF-1/WFV, Terra/MODIS, Himawari-8/AHI AODs derived from different retrieval algorithms. Figures 9 and 10 depict the spatial distributions of AOD_WMRT (500 m), AOD_WORT (500 m), AOD_MODT (10 km), AOD_MODB (10 km), and AOD_AHI (0.05 deg) on 12 July 2017 and 1 January 2018, respectively. The heavily polluted areas were mainly concentrated on the southeast part of BTH region, which were obviously different from the northwest mountainous areas. The spatial distribution of our AOD_WORT retrievals is generally consistent with those of other four satellite AODs. The high values of AOD (>0.8) are mainly distributed in the southeast of BTH region, and zones with the low AODs are located in the northwest parts. In the northwestern BTH region with lower AOD values area, AOD_WORT, and AOD_MODB are much smoother than other three AOD results. It is notable that the invalid value ratio of GF-1/WFV AODs is much lower than those of Terra/MODIS, Himawari-8/AHI ones, especially for AOD_MODT and AOD_AHI, which may be caused by the difference in spatial resolution, overpassing time, and retrieval algorithm. Such a situation of a high invalid value ratio was particularly serious for the AOD_MODT over the southeastern BHT region with air pollution, which made the status over this region cannot be captured for the analysis. And some pixels covered with snow and ice in Figure 9b,c and Figure 10b,c were also set as invalid values the northwest mountains. However, for the retrieved GF-1/WFV AOD with a spatial resolution of 500 m, both AOD_WORT and AOD_WMRT can provides more abundant details than MODIS and AHI in the spatial distribution, and can more accurately reflect the AOD spatial variation over the polluted region. serious for the AOD_MODT over the southeastern BHT region with air pollution, which made the status over this region cannot be captured for the analysis. And some pixels covered with snow and ice in Figures 9b,c and 10b,c were also set as invalid values the northwest mountains. However, for the retrieved GF-1/WFV AOD with a spatial resolution of 500 m, both AOD_WORT and AOD_WMRT can provides more abundant details than MODIS and AHI in the spatial distribution, and can more accurately reflect the AOD spatial variation over the polluted region.  The reflectance in bright areas (mainly include urban, desert, bare land, and arid/semiarid areas with sparse or little vegetation coverage) is more difficult to estimate accurately [14,15]. Central Beijing is a typically bright urban area with a dense population and extensive industrial activity. As shown in Figures 9 and 10, the AOD distribution in Beijing is weakly affected by the bright surface. The bright surface of urban area does not cause high value distribution of retrieval results. The spatial distribution of AOD in urban and non-urban areas is relatively smooth, and the removal of surface information is more accurate. So, it is indicated that the high AODs were not caused by the bright surface of urban area in our retrieval algorithm.
Because the LSR of the current month cannot be obtained for the near-real-time AOD retrieval in the operational application, in our study, the land surface reflectance constructed by the ORT method was a monthly dataset transformed from the MODIS LSR data during the period from 2016 to 2018. In practice, the monthly LSR based on the MODIS LSRs of last three-year will be used to retrieve the near-real-time GF-1/WFV AOD of the current year, where the month of GF-1/WFV data is the same as that of monthly MODIS LSR. And the LSRs of the current month will be continuously updated into our monthly LSR dataset. According to the true color image on 28 September 2019 (Figure 11a), there was a serious pollution event in the southeast of the BTH region, which can be directly and significantly reflected in our retrieved GF-1/WFV AOD_ WORT result (Figure 11b). The GF-1/WFV AOD_WORT can also be used to strictly distinguish the heavy polluted regions on an urban-scale, which is more useful for the regional control and governance. The reflectance in bright areas (mainly include urban, desert, bare land, and arid/semi-arid areas with sparse or little vegetation coverage) is more difficult to estimate accurately [14,15]. Central Beijing is a typically bright urban area with a dense population and extensive industrial activity. As shown in Figures 9 and 10, the AOD distribution in Beijing is weakly affected by the bright surface. The bright surface of urban area does not cause high value distribution of retrieval results. The spatial distribution of AOD in urban and non-urban areas is relatively smooth, and the removal of surface information is more accurate. So, it is indicated that the high AODs were not caused by the bright surface of urban area in our retrieval algorithm. Because the LSR of the current month cannot be obtained for the near-real-time AOD retrieval in the operational application, in our study, the land surface reflectance constructed by the ORT method was a monthly dataset transformed from the MODIS LSR data during the period from 2016 to 2018. In practice, the monthly LSR based on the MODIS LSRs of last three-year will be used to retrieve the near-real-time GF-1/WFV AOD of the current year, where the month of GF-1/WFV data is the same as that of monthly MODIS LSR. And the LSRs of the current month will be continuously updated into our monthly LSR dataset. According to the true color image on 28 September 2019 ( Figure  11a), there was a serious pollution event in the southeast of the BTH region, which can be directly and significantly reflected in our retrieved GF-1/WFV AOD_ WORT result (Figure 11b). The GF-1/WFV AOD_WORT can also be used to strictly distinguish the heavy polluted regions on an urban-scale, which is more useful for the regional control and governance.

AOD Validation Using AERONET Data
Ground-based AOD measurements of five AERONET sites (i.e., Beijing, XiangHe, Beijing-CAMS, Beijing-PKU, and Beijing-RADI sites) were extracted for the AOD validation. The error statistics, including the root mean square error (RMSE) and the mean absolute error (MAE), were verified by comparing the satellite-retrieved AODs (AOD_WORT and AOD_WMRT) with the AERONET measurements (AOD_AERO). RMSE was selected to measure the differences between AOD retrievals and ground-based measurements, which was sensitive to both systematic and random errors. MAE could better reflect the actual situation of AOD retrievals error. RMSE and MAE can be written as follows, respectively:

AOD Validation Using AERONET Data
Ground-based AOD measurements of five AERONET sites (i.e., Beijing, XiangHe, Beijing-CAMS, Beijing-PKU, and Beijing-RADI sites) were extracted for the AOD validation. The error statistics, including the root mean square error (RMSE) and the mean absolute error (MAE), were verified by comparing the satellite-retrieved AODs (AOD_WORT and AOD_WMRT) with the AERONET measurements (AOD_AERO). RMSE was selected to measure the differences between AOD retrievals and ground-based measurements, which was sensitive to both systematic and random errors. MAE could better reflect the actual situation of AOD retrievals error. RMSE and MAE can be written as follows, respectively: Both of AOD_WMRT and AOD_WORT retrievals were matched with AERONET AOD measurements following the rule mentioned in Section 2.2. The validation results of AOD_WMRT and AOD_WORT from 2016 to 2018 were shown in Figure 12. Compared with AOD_WMRT (R = 0.79), AOD_WORT with a higher correlation coefficient (R = 0.88) showed a better relationship with ground measurements. AOD_WMRT was significantly overestimated for low AOD values, which might be due to the undervaluation of LSR. However, our AOD_ WORT retrieval algorithm significantly modified such overestimation of AOD_WMRT. Especially, when the AOD is smaller than 0.2, a better agreement between satellite AOD retrievals and AERONET observations for AOD_ WORT than that for AOD_WMRT. Separately, the AOD_WORT retrievals on Beijing, Beijing-CAMS, and XiangHe sites were much better agreement with AERONET observations. showed a better relationship with ground measurements. AOD_WMRT was significantly overestimated for low AOD values, which might be due to the undervaluation of LSR. However, our AOD_ WORT retrieval algorithm significantly modified such overestimation of AOD_WMRT. Especially, when the AOD is smaller than 0.2, a better agreement between satellite AOD retrievals and AERONET observations for AOD_ WORT than that for AOD_WMRT. Separately, the AOD_WORT retrievals on Beijing, Beijing-CAMS, and XiangHe sites were much better agreement with AERONET observations. To testify whether the LSR dataset based on 2016-2018 MODIS LSR data is proper for retrieving the GF-1/WFV AOD of 2019, the same validation work for 2019 AOD_WORT (Figure 13) was done as for 2016-2018 AOD_WORT ( Figure 12). Figures 13  and 14 indicate a better performance for our proposed algorithm in this paper. R of 2019 AOD_WORT reaches up to 0.94, especially R of the XiangHe site is the largest one among that of all sites, as high as 0.98. Both of RMSE and MAE for AOD_WORT are much lower than those for AOD_WMRT. For R, the largest one is 2019 AOD_WORT, and following is 2016-2018 AOD_WORT, and the smallest one is 2016-2018 AOD_WMRT. On the contrary, both of RMSE and MAE for 2019 AOD_WORT are smaller than those of 2016-2018 AOD_WORT and AOD_WMRT. Overall, for all 5 AERONET sites over BTH region, a relatively high R (R > 0.88) and low RMSE and MAE (RMSE < 0.13, MAE < 0.11) appear between AOD_ WORT and AOD_AERO.
Overall, although there were random uncertainties in the dataset building of GF-1/WFV surface reflectance led by the MODIS LSR errors, the proposed ORT method could reduce such uncertainties by dealing strictly with MODIS surface reflectance and taking To testify whether the LSR dataset based on 2016-2018 MODIS LSR data is proper for retrieving the GF-1/WFV AOD of 2019, the same validation work for 2019 AOD_WORT (Figure 13) was done as for 2016-2018 AOD_WORT ( Figure 12). Figures 13 and 14 indicate a better performance for our proposed algorithm in this paper. R of 2019 AOD_WORT reaches up to 0.94, especially R of the XiangHe site is the largest one among that of all sites, as high as 0.98. Both of RMSE and MAE for AOD_WORT are much lower than those for AOD_WMRT. For R, the largest one is 2019 AOD_WORT, and following is 2016-2018 AOD_WORT, and the smallest one is 2016-2018 AOD_WMRT. On the contrary, both of RMSE and MAE for 2019 AOD_WORT are smaller than those of 2016-2018 AOD_WORT and AOD_WMRT. Overall, for all 5 AERONET sites over BTH region, a relatively high R (R > 0.88) and low RMSE and MAE (RMSE < 0.13, MAE < 0.11) appear between AOD_ WORT and AOD_AERO.
Overall, although there were random uncertainties in the dataset building of GF-1/WFV surface reflectance led by the MODIS LSR errors, the proposed ORT method could reduce such uncertainties by dealing strictly with MODIS surface reflectance and taking the averaged LSR values. Our GF-1/WFV AOD retrieval algorithm performed well in the improvement of the retrieval results, and the AOD inversions showed a good correlation and accuracy by comparing with multi-source AOD products. Especially, the AOD_WORT is superior to AOD_WMRT for low AOD values. Whether the LSR dataset built by 2016-2018 MODIS LSR data was used to retrieve the GF-1/WFV AOD of 2016-2018 or 2019, both the AOD retrievals of the two periods were shown with high accuracy.

Conclusions
In this study, a new DB AOD retrieval algorithm was developed for GF-1/WFV, involving a cloud screening method based on the only four bands of GF-1/WFV and an optimal reflectivity technique method for building a monthly LSR dataset transmitted from MODIS LSR. GF-1/WFV AOD_WORT from 2016 to 2019 was retrieved and validated by multi-source satellite AOD products and AERONET ground-based measurements. Detailed validation and comparison results were shown as follows: 1. In terms of spatial distribution, our AOD_WORT results were roughly consistent with Terra/MODIS and Himawari-8/AHI AOD retrievals based on different AOD re-

Conclusions
In this study, a new DB AOD retrieval algorithm was developed for GF-1/WFV, involving a cloud screening method based on the only four bands of GF-1/WFV and an optimal reflectivity technique method for building a monthly LSR dataset transmitted from MODIS LSR. GF-1/WFV AOD_WORT from 2016 to 2019 was retrieved and validated by multi-source satellite AOD products and AERONET ground-based measurements. Detailed validation and comparison results were shown as follows: 1. In terms of spatial distribution, our AOD_WORT results were roughly consistent with Terra/MODIS and Himawari-8/AHI AOD retrievals based on different AOD re-

Conclusions
In this study, a new DB AOD retrieval algorithm was developed for GF-1/WFV, involving a cloud screening method based on the only four bands of GF-1/WFV and an optimal reflectivity technique method for building a monthly LSR dataset transmitted from MODIS LSR. GF-1/WFV AOD_WORT from 2016 to 2019 was retrieved and validated by multi-source satellite AOD products and AERONET ground-based measurements. Detailed validation and comparison results were shown as follows: 1.
In terms of spatial distribution, our AOD_WORT results were roughly consistent with Terra/MODIS and Himawari-8/AHI AOD retrievals based on different AOD retrieval algorithms. GF-1/WFV AOD with a higher spatial resolution of 500 m can provide more abundant details of spatial variation than Terra/MODIS and Himawari-8/AHI.

3.
The near-real-time inversion of AOD based on the LSR dataset of last three years also performed excellently both in terms of the spatial distribution and the quantitative comparison with AERONET observations. The high correlation with AERONET AOD (R ORT = 0.94), the low RMSE (RMSE ORT = 0.11) and the low MAE (MAE ORT = 0.09) could certificate the reliability of our AOD_WORT retrieval algorithm.

4.
Although there was no detailed comparison of inversion effects between rural and urban backgrounds, the smoothing effect on the spatial distribution could confirm that the presence of urban bright surfaces would not cause high values or noises.
Our GF-1/WFV AOD inversion study can lay a foundation for other similar sensors aboard Gaofen satellite series (such as GF-1B, GF-1C, GF-2, and GF-6), in terms of the solution of the inversion problem, the optimization of the inversion algorithm, and the processing and comparison of the inversion verification.