An Improved Aerosol Optical Depth Retrieval Algorithm for Multiangle Directional Polarimetric Camera (DPC)

: The DPC is a multiangle sensor that detects atmospheric parameters. However, the retrieval of high-precision and high-spatial-resolution aerosol products from the DPC remains a great challenge due to the ill-posed nature of the problem. Thus, a novel aerosol optical depth (AOD) retrieval algorithm was proposed using visible surface reﬂectance relationships (VISRRs). The VISRR algorithm accounts for the surface anisotropy and needs neither a shortwave infrared band nor a surface reﬂectance database that can retrieve AOD over dark and bright land cover. Firstly, moderate-resolution imaging spectroradiometer (MODIS) surface reﬂectance (MYD09) products were used to derive the preceding surface reﬂectance relationships (SRRs), which are related to surface types, scattering angle, and normalized difference vegetation index (NDVI). Furthermore, to solve the problem of the NDVI being susceptible to the atmosphere, an innovative method based on an iterative atmospheric correction was proposed to provide a realistic NDVI. The VISRR algorithm was then applied to the thirteen months of DPC multiangle data over the China region. AOD product comparison between the DPC and MODIS showed that they had similar spatial distribution, but the DPC had both high spatial resolution and coverage. The validation between the ground-based sites and the retrieval results showed that the DPC AOD performed best, with a Pearson correlation coefﬁcient (R) of 0.88, a root mean square error (RMSE) of 0.17, and a good fraction (Gfrac) of 62.7%. Then, the uncertainties regarding the AOD products were discussed for future improvements. Our results revealed that the VISRR algorithm is an effective method for retrieving reliable, simultaneously high-spatial-resolution and full-surface-coverage AOD data with good accuracy.


Introduction
Aerosol particles are widely suspended in the Earth's atmosphere, which greatly impacts the Earth's radiation balance and climate change by scattering and absorbing solar light. Meanwhile, aerosol particles can cause adverse effects on the atmospheric environment, human health, satellite image quality, etc. Thus, monitoring and determining aerosols' optical, physical, and chemical properties are necessary for different research fields. There are two main aerosol monitoring methods based on the principle of aerosols' extinction effect on solar radiation. One uses ground-based sun photometers, and the other uses satellite-based sensors. Aerosol products are often measured via ground-based observations with high precision that are also used to validate satellite-based aerosol and surface contrast is sufficient in a certain size window. Thus, it requires a reduction in spatial resolution, and it is hard to retrieve the aerosol properties at the intrinsic spatial resolution [29]. For polarized sensor aerosol retrieval, the surface polarized reflectance (SPR) estimation is usually derived using the bidirectional polarization distribution function (BPDF) model [14]. However, the BPDF may yield relatively large estimation errors in some surface types. To solve this problem, Ge et al. [30] provided an alternative approach that utilizes the spectral neutrality characteristic of SPR as a constraint to retrieve higherprecision fine-mode aerosol products. In addition, a statistically optimized algorithm called generalized retrieval of aerosol and surface properties (GRASP) was designed for characterizing atmosphere and surface properties simultaneously from a diversity of remote sensing measurements [31]. The aerosol products of GRASP have been validated as highquality at a global and regional scale [32,33].
Regarding prior knowledge of the aerosol components, the common method is to construct typical aerosol models using ground-based observations. Aerosol model properties usually include the imaginary and real part of the refractive index, parameters of size distribution, and the single-scattering albedo. For example, Dubovik et al. [34] derived the aerosol absorption and other properties in several key locations based on the long time-series data of global AERONET sites. Omar et al. [35] derived six aerosol models, i.e., biomass burning, urban industrial pollution, rural background, polluted marine, dirty pollution, and desert dust; Levy et al. [36] used the cluster analysis method and derived four aerosol models that varied according to the location and season, i.e., moderately absorbing, absorbing, non-absorbing, and spheroid; and Lee et al. [37] derived six aerosol models for East Asia and showed that the models were strongly affected by their sources. Meanwhile, these aerosol models have been widely used in satellite remote sensing aerosol retrieval [7,19,23,38,39]. However, due to the sites being sparse and distributed unevenly, the aerosol models based on AERONET may not be representative of China. Thus, Li et al. [40] used the long time-series of SONET observations to derive fundamental aerosol models for China.
When the two basic problems have been solved, a look-up table (LUT) method is often used in operational retrieval algorithms, such as the MODIS, VIIRS, AATSR, and MISR, to save computation time. Meanwhile, aerosol products retrieved from different sensors have been well-validated globally [10,41,42]. However, these aerosol products may have the shortcomings of low spatial resolution (e.g., the aerosol products (MXD04) of DT and DB have a spatial resolution of 10 km) and low precision for local and regional research [43][44][45][46][47][48].
In May 2018, the multiangle multispectral polarized sensor directional polarimetric camera (DPC) onboard the Chinese GaoFen-5 (GF-5 01) satellite was successfully launched [13]. The DPC is equipped with three polarized bands (490, 670, and 865 nm) and five non-polarized bands (443, 565, 763, 765, and 910 nm). The spatial resolution and swath width of the DPC are 3.3 and 1850 km, respectively, and the maximum observation angle of the DPC can be up to 12. The DPC is evidently a powerful tool for monitoring atmospheric parameters such as aerosols, clouds, and water vapor. Unfortunately, due to the failure of the solar panel, DPC observation only lasted until April 2020. It provided us with 13 months of continuous observational data, that is, from March 2019 to March 2020. However, as described above, the current operational algorithms have some shortcomings. Furthermore, most of them do not consider surface anisotropy. Thus, to retrieve highspatial-resolution and high-precision aerosol products over both dark and bright land cover based on multiangle intensity data from the DPC, a novel aerosol retrieval algorithm based on the visible spectral surface reflectance relationships (VISRR) method was proposed (the detailed VISRR algorithm is introduced in Section 3.3). In addition to the DPC/GF-5 (01), the VISRR algorithm can also be applied to similar sensors, such as the advanced DPCs onboard the Chinese GF-5 (02), carbon monitoring (CM-1), and atmospheric environmental monitoring (DQ-1 and DQ-2) satellites and the European multi-view multi-channel multi-polarization imaging (3MI) onboard EPS-SG, etc. [49]. When combined with the fine-model AOD, e.g., our previous fine-model aerosol algorithm 'spectral neutrality of Remote Sens. 2022, 14, 4045 4 of 25 surface polarized reflectance' (SNOSPR), based on multiangle polarized data [30], the fine-model aerosol fraction (FMF) can be calculated, which can characterize the weight of aerosols from anthropogenic and natural sources. This paper is structured as follows: In Section 2, the study area and datasets are presented. Section 3 introduces the basic principles of aerosol retrieval based on satellite remote sensing and the VISRR algorithm. Then, in Section 4, the spatial distribution of daily and average DPC AOD products is presented and compared with that of MODIS products. Then, all aerosol products are validated with ground-based (i.e., AERONET and SONET) products, and the uncertainties regarding the AOD are analyzed and discussed. Section 5 presents the conclusions.

Study Area
The territory of China was the study area, and its land cover types are shown in Figure 1. According to the MODIS Land Cover Type (MCD12Q1) International Geosphere-Biosphere Programme (IGBP) classification, China has 17 land cover types (as shown in Table 1). The spatial resolution was 500 m [50]. It can be seen that China comprises all surface types, with vast deserts and grasslands distributed in the northwest, croplands mainly located in the east and northeast, forests in the south and northeast, and urban cities mainly distributed in the east and along the coastal area. Furthermore, China is currently the biggest developing country globally and has the largest population; with its rapid development in the past few decades, China has become one of the primary sources of aerosols globally. Thus, satellite aerosol remote sensing over China is necessary and has always been a research hotspot. However, due to the complexity of the aerosol sources and land cover, aerosol retrieval over China is still a great challenge.

Elevation Data
Global Multiresolution Terrain Elevation Data 2010 (GMTED2010) is a global elevation model developed by the U.S. Geological Survey (USGS, https://usgs.gov/) (accessed on 15 August 2022) and National Geospatial Intelligence Agency (NGA), which has three spatial resolutions of about 1 km, 500 m, and 250 m [52]. In this study, GMTED2010 data were used to determine the altitude and correct the influence of Rayleigh scattering for the DPC.

Spectral Library
The USGS spectral library is a reference database containing seven categories of reflectance spectra (artificial materials, coatings, liquids, minerals, organic compounds, soils and mixtures, and vegetation) measured using laboratory, field, and airborne spectrometers [53]. The Spectral Library Version07b (Splib07b) datasets were used to analyze the characteristic variety of the surface spectra between different bands, in order to adjust the spectral response function (SRF) differences between DPC and MODIS for corresponding bands.  [1], and the uncertainty of the VIS and NIR AERONET AOD is estimated to be 0.01 for Version 3 (V3) products. However, AERONET sites are rare and distributed unevenly over China. Thus, the Chinese ground-based aerosol observation network SONET (http://www.sonet.ac.cn/en/index.php) (accessed on 15 August 2022) was built by the Chinese Academy of Sciences in collaboration with universities and institutes. The sites of SONET are widely distributed in the typical area of China [2]. In this study, twenty-two AERONET and SONET sites were chosen to validate the DPC and MODIS aerosol products. The spatial distribution and detailed information related to these sites are presented in Figure 1 and Table 2. In addition, these sites were divided into five surface types according to MCD12Q1 [54].

Surface BRDF Model
This study adopted a linear semi-empirical kernel-driven RossThick LiSparseReciprocal (RTLSR) BRDF model to describe the anisotropic SR [55,56]. The RTLSR model has been widely used in many types of research and has obtained good results [57,58]. Meanwhile, the RTLSR model was adopted to produce MODIS and VIIRS global operational BRDF products [59,60]. The equation of the RTLSR model is as follows: where λ is the wavelength; θ s , θ v , and ∅ are solar zenith angle, sensor zenith angle, and relative azimuth angle, respectively; R SR represents the surface directional reflectance; K iso , K vol , and K geo represent the isotropic, volumetric, and geometric kernels, respectively, where K iso = 1, and the expression of K vol and K geo , which are fixed functions and only dependent on the solar and viewing geometry, can be found in Wanner et al. [61]; and f iso (λ), f geo (λ), and f vol (λ) are the coefficients of the three kernels, which were varied according to the wavelength and surface characteristics. Furthermore, Equation (1) could be separated into two terms: where ρ(λ) = f iso (λ) indicates the "reflectance magnitude", which is governed by surface type microphysical properties and changes rapidly with wavelength and time; α 1 and α 2 represent the geometric factor and volumetric factor, respectively; α 1 (λ) = f vol (λ) f iso (λ) ; and α 2 (λ) = f geo (λ) f iso (λ) . The part in the square bracket is the BRDF shape (BRDF s ) function (Equation (3)), which is based on the macroscopic structure of the surface types, remains nearly independent of the wavelength, and changes slowly in a short time [62]:

Atmospheric Radiative Transfer Model
Overlying the Lambertian homogeneous surface and under cloud-free, plane-parallel atmosphere conditions, the satellite sensor received reflectance (R TOA ) at the top of atmosphere (TOA) can be written as: where τ is the optical depth; T g indicates the gaseous transmission caused by gases such as column water vapor (CWV) and ozone (O 3 ); R atm refers to the aerosol and molecular intrinsic reflectance; T ↓ and T ↑ are transmission; S is the spherical albedo of the atmosphere; and R SR is the angular SR. SR is an anisotropic Lambertian surface hypothesis can affect the accuracy of the R TOA simulation, especially for multiangle data, ultimately transferring the error to the aerosol retrieval [63]. To avoid this problem in the present study, a fast, accurate, non-Lambertian atmospheric radiative transfer function considering the anisotropy of the surface based on the four-stream theory was adopted. This function has high precision, with mean relative differences in the spectral (range from UV to NIR) TOA simulation of less than 0.7% for different surface types [64]. This forward model has been widely used for aerosol and surface parameter retrieval [65][66][67], and the equation can be expressed as: ; t s (θ s/v ) refers to the diffuse transmission; R is the reflectance matrix; and |R| is the determinant of R.
where R DHR represents the directional-hemispherical reflectance (DHR), i.e., black-sky albedo (R bs ), which describes the diffuse reflection of the incoming direct beam over the hemisphere; R HDR refers to the hemispherical-directional reflectance (HDR), which describes the direct reflection of the incoming diffuse radiation from the whole hemisphere; and R BHR is the bi-hemispherical reflectance (BHR), i.e., white-sky albedo. When the coefficients of the three kernels are obtained, the SR for any angle can be calculated with Equation (1). Meanwhile, the R DHR and R BHR can be computed using the following Equations (8)-(11) [67]: where H k and g ik (i = 0, 1, 2, 3; k = iso, vol, geo) are the regression coefficients listed in Table 3. Table 3. Regression coefficients used to calculate black-sky and white-sky albedo.

VISRR Algorithm
VISRR is a novel AOD retrieval algorithm, which uses the SRRs between VIS bands as constraints and considers the influence of surface anisotropy. VISRR can simultaneously retrieve AOD over dark and bright land covers; the algorithm derivation and retrieval strategies are described in detail below.

SR Characteristics Analysis
Firstly, we used the Splib07b spectral library to analyze which bands were suitable for aerosol retrieval. Figure 2 shows the SRFs of the DPC in the 443, 490, 565, 670, and 865 nm bands (red dotted line) and the corresponding MODIS b9, b10, b5, b1, and b2 bands (gray dotted line). Curves of different colors represent the SR spectral curves of typical surface materials (i.e., blackbrush, sand, brick, and acmite) from Splib07b. From the picture, it can be seen that the surface contribution of VIS is much lower than that of NIR; the SR of non-green vegetation (i.e., sand and brick) is higher than that of green vegetation (i.e., blackbrush); and for green vegetation, there exist two high-reflection peaks near 565 and 865 nm. These features mean that VIS bands are more suitable for aerosol retrieval than NIR bands; aerosol retrieval over green vegetation is easier than over non-green vegetation; and the green band (565 nm) may not be suitable for aerosol retrieval over green vegetation cover. Furthermore, we analyzed the SRRs between the VIS and NIR bands; Figure 3a,b show the SRRs between 490 and 670 nm and between 490 and 865 nm, respectively. It can be seen that the SRRs between VIS bands are more stable than those between VIS and NIR bands, which is consistent with Levy's research [19]. This phenomenon can be explained by the fact that most ground objects have similar spectral reflectance changes in the VIS spectrum range but large changes in the NIR or SWIR bands. Thus, the 443, 490, and 670 nm bands were chosen to build the SRRs in this study.
in the VIS spectrum range but large changes in the NIR or SWIR bands. Thus, the and 670 nm bands were chosen to build the SRRs in this study.

SR Relationship Analysis
Within a satellite sensor's instantaneous field of view, individual pixels det the scanning device often include more than one object, especially those with a resolution at the kilometer level. Thus, the determination of SRRs requires real data. Two common methods are used to construct the target sensor's prior know spectral SR for aerosol retrieval. One is to apply a minimum reflectivity techniqu or atmospheric correction to the target sensor's data [7,19]. This method is suit target sensors with long-time-series data that can provide enough samples. Th method is to use the mature SR products from other satellite sensors, such as [69,70]. This method is generally applicable to new sensors. In this study, owin limitation of the relatively small amount of DPC data, the second method was cho three years (2016-2018) of MYD09 SR products pertaining to China and its surr The solid lines show the spectral reflectance of typical surface materials in Splib07b, i.e., blackbrush, sand, brick, and acmite, representing vegetation, soil, artificial material, and minerals, respectively. The abscissa represents the wavelength, the left ordinate represents the spectral response function value, and the right represents the surface reflectance value.
Remote Sens. 2022, 14, x FOR PEER REVIEW 9 of 2 in the VIS spectrum range but large changes in the NIR or SWIR bands. Thus, the 443, 490 and 670 nm bands were chosen to build the SRRs in this study.

SR Relationship Analysis
Within a satellite sensor's instantaneous field of view, individual pixels detected by the scanning device often include more than one object, especially those with a spatia resolution at the kilometer level. Thus, the determination of SRRs requires real satellite data. Two common methods are used to construct the target sensor's prior knowledge o spectral SR for aerosol retrieval. One is to apply a minimum reflectivity technique [17,68 or atmospheric correction to the target sensor's data [7,19]. This method is suitable fo target sensors with long-time-series data that can provide enough samples. The othe method is to use the mature SR products from other satellite sensors, such as MODIS [69,70]. This method is generally applicable to new sensors. In this study, owing to th limitation of the relatively small amount of DPC data, the second method was chosen, and three years (2016-2018) of MYD09 SR products pertaining to China and its surrounding areas ( Figure 1) were used to derive the DPC equivalent spectrum SR over different sur face types. To ensure the accuracy of the SR, only the pixels with an AOD at 470 nm

SR Relationship Analysis
Within a satellite sensor's instantaneous field of view, individual pixels detected by the scanning device often include more than one object, especially those with a spatial resolution at the kilometer level. Thus, the determination of SRRs requires real satellite data. Two common methods are used to construct the target sensor's prior knowledge of spectral SR for aerosol retrieval. One is to apply a minimum reflectivity technique [17,68] or atmospheric correction to the target sensor's data [7,19]. This method is suitable for target sensors with long-time-series data that can provide enough samples. The other method is to use the mature SR products from other satellite sensors, such as MODIS [69,70]. This method is generally applicable to new sensors. In this study, owing to the limitation of the relatively small amount of DPC data, the second method was chosen, and three years (2016-2018) of MYD09 SR products pertaining to China and its surrounding areas ( Figure 1) were used to derive the DPC equivalent spectrum SR over different surface types. To ensure the accuracy of the SR, only the pixels with an AOD at 470 nm smaller than 0.2 and the highest quality assurance (QA) of MYD09 were chosen for analysis.
However, two problems needed to be solved beforehand, since DPC and MODIS differ in terms of spatial resolution (3.3 versus 1 km) and SRFs (as shown in Figure 2). To make the spatial resolution consistent, about 3 × 3 pixels of MYD09 and MYD03 were averaged at 3.3 km. As shown in Figure 2, the SRFs of DPC and MODIS do not exactly match; if we did not consider these differences, it would have led to errors in the derived equivalent spectral SR of DPC and affected the accuracy of the aerosol retrieval [70]. We adopted the singular-value decomposition (SVD) method to predict the DPC equivalent spectral SR from the corresponding MODIS channels to solve this problem. The steps were as follows: (1) we used the SRFs and the spectral curves of typical surface materials from Splib07b to calculate the DPC and MODIS spectral SR; (2) we calculated the SVDs, which included the singular vectors V modis and V DPC , based on the calculated spectral SR from step (1); (3) we used the V modis and MYD09 SR products to calculate the conversion coefficient vector C; and (4) the DPC equivalent spectrum SR could be calculated using the C and the V DPC . Detailed information on this method can be found in Sayer et al. [70].
After obtaining the equivalent spectrum SR and geometric information on the DPC over different surface types (except for permanent snow, ice, and water), we analyzed the SRRs between 443, 490, and 670 nm. The normalized difference vegetation index (NDVI) is an important adjustment parameter to estimate the SR for aerosol retrieval [19,20,71], and it can be calculated by the SR of 670 and 865 nm using Equation (12): The statistical results of the SRRs are shown in Figure 4. It can be seen that the Pearson correlation coefficient (R) of 443 versus 670 nm ranged from 0.82 to 0.98, the root mean square error (RMSE) ranged from 0.003 to 0.012, and the average SR ratios (K 443_670 ) over different surface types varied from 0.43 to 0.63. For 490 and 670 nm, the R ranged from 0.90 to 0.99, the RMSE ranged from 0.002 to 0.011, and the K 490_670 over different surface types varied from 0.51 to 0.68. These statistical parameters showed that: (1) the SRRs between VIS bands had a high correlation, which was consistent with the conclusion of Section 3.3.1; (2) the SRRs between VIS bands needed to be divided according to the surface types due the substantial differences; (3) the correlation of SRRs between 490 and 670 nm was higher than that between 443 and 670 nm overall; and (4) the NDVI could indicate the surface type and the change in the SR over time to a certain extent-for example, the NDVI of evergreen vegetation (e.g., Figure 4(1,2)) usually has a high value, while over bare soil (e.g., Figure 4 (16)) it has a low value. When the surface types vary significantly with the seasons, the NDVI also has a wider range of changes (e.g., Figure 4(10,13)). It is necessary to consider surface types and NDVI when using the SRRs between VIS bands as constraints for aerosol retrieval.
In addition, to make the SRRs between VIS bands more accurate, the SCA was also taken into account when determining the K 443_670 and K 490_670 over different surface types. The final values of K 443_670 and K 490_670 over each typical surface type are shown in Tables 4-8 (mixed forests, grasslands, croplands, urban and built-up lands, and barren land). The NDVI was set from 0 to 1 with an interval of 0.2; the SCA was set from 60 • to 100 • and from 100 • to 180 • with an interval of 20 • . The ratios over other surface types (i.e., (1)-(4), (6)-(9), (14), and (18)- (21)) are shown in Supplemental Materials Section S2, Tables S1-S10. In addition, to make the SRRs between VIS bands more accurate, the SCA was also taken into account when determining the 443_670 and 490_670 over different surface types. The final values of 443_670 and 490_670 over each typical surface type are shown in Tables 4-8 (mixed forests, grasslands, croplands, urban and built-up lands, and barren land). The NDVI was set from 0 to 1 with an interval of 0.2; the SCA was set from 60° to 100° and from 100° to 180° with an interval of 20°. The ratios over other surface types (i.e., (1)-(4), (6)-(9), (14), and (18)-(21)) are shown in Supplemental Materials Section S2, Tables S1-S10.

Aerosol Models
The aerosol model is another essential prior knowledge component, and its accuracy plays a vital role in the satellite remote sensing of aerosols [72]. Several global and regional aerosol models derived from AERONET measurements [35][36][37] have been used for many research satellite aerosol retrieval methods. However, these models may not be very representative in China, where ground-based sites are sparse and distributed unevenly. For example, the aerosol models used in the VIIRS enterprising processing system (EPS) aerosol retrieval algorithm have been validated as unsatisfactory in China [44]. Thus, based on the long-time-series data of the SONET, Li et al. [40] used the K-means cluster analysis approach and obtained ten fundamental aerosol models for China, including five typical fine-particle aerosol models ((1) urban polluted, (2) secondary polluted, (3) combined polluted, (4) polluted fly ash, and (5) continental background) and five coarse models ((6) summer fly ash, (7) winter fly ash, (8) primary dust, (9) transported dust, and (10) background dust). Li et al. [40] also provided the probability of the appearance of fine-mode aerosols and coarse-mode aerosols. In this study, we adopted the six combinations with the highest probability of appearance for aerosol retrieval, i.e., urban polluted and summer fly ash, secondary polluted and summer fly ash, combined polluted and winter fly ash, continental background and background dust, continental background and winter fly ash, and continental background and transported dust. The microphysical parameters of these aerosol models can be found in Li et al. [40].

Retrieval Scheme
In this study, the VISRR algorithm was implemented based on the precalculated LUTs, which were built by the second simulation of a satellite signal in the solar spectrum vector (6SV) radiative transfer (RT) model [73]. The input parameters of the 6SV RT model included the SRFs of DPC 443, 490, 670, and 865 nm; six combinations of fine-mode aerosols and coarse-mode aerosols, as shown in Section 3.3.3; θ s and θ v ranging from 0 • to 80 • with an interval of 6 • ; ∅ ranging from 0 • to 180 • with an interval of 12 • ; 11 AOD values (550 nm)-0, 0.01, 0.25, 0.5, 1, 1.5, 2, 2.5, 3, 4, 5; and elevation set to sea level. Then, atmospheric parameters (i.e., R atm , → T (θ s ), → T (θ v ), t s (θ s ), and t s (θ v ) in Equation (4)) for various observational and atmospheric scenarios were built in the LUTs. Similarly, the LUTs of CWV and O 3 gas transmission were also built by the 6SV RT model. The value of CWV was set in the range of 0 to 7 g/cm 2 with an interval of 0.5, and the value of O 3 was set in the range of 0 to 5 cm-atm with an interval of 0.5. In this study, the prior knowledge of BRDF shape (Equation (3)) was assumed to be unchanged for one month and independent of bands [62,65]. When the LUTs and related auxiliary data were ready, the VISRR algorithm could be performed based on the multiangle multispectral DPC intensity data. All procedures were applied to the individual pixels, and the output result at each pixel was recorded. The schematic flowchart is shown in Figure 5, and the specific retrieval steps are described below.
The first step was the identification of the land surface type for each pixel using the MCD12Q1 products; second was the recognition of cloud-contaminated pixels, as only cloud-free and land (except for permanent snow, ice, or inland water) pixels were further processed in the next step; and third was to match and extract the auxiliary data according to the longitude, latitude, the corresponding geometric factor α 1 and volumetric factor α 2 of the BRDF shape, the daily O 3 and CWV, and the elevation from the month-averaged MCD43D19-20, MYD08_D3, and GMTED2010 products. The global mean values were adopted if the O 3 and CWV were missing from MYD08_D3, as shown in Levy et al. [74]. (2) We utilized the Rayleigh scattering correction method proposed by Fraser et al. [75] to remove its contribution to the reflectance at TOA [19] and then obtained the gas corrected reflectance (R Meas_gasc TOA (λ)).
(3) Considering the non-Lambertian characteristics of the surface, the atmospheric correction (AC) was implemented to calculate accurate directional surface reflectance R SR (λ). The AC equation is Equation (15), and the detailed derivation process is shown in Supplementary Section S3.
where the definitions of F 1 , F 2 , F 3 , F 4 can be found in Supplementary Section S3. Then, the three kernel coefficients f iso (λ), f geo (λ), and f vol (λ) of the RTLSR BRDF model and the directional surface reflectance R SR (λ) could be obtained as follows: (c) Determination of the initial NDVI value and prior knowledge of K 443_670 and K 490_670 The NDVI is an important adjustable parameter for estimating SR and has been widely used for aerosol retrieval. However, the NDVI usually cannot be obtained before aerosol retrieval due to the effects of the atmosphere. Thus, some researchers directly use the reflectance at TOA to calculate the NDVI (NDVI_TOA) [76,77] or use the Rayleigh corrected reflectance at TOA to calculate the NDVI (NDVI_RC) to reduce the impact of the atmosphere [20,78]. However, NDVI_TOA and NDVI_RC are still susceptible to aerosols and rapidly decrease with increasing AOD [79][80][81]. The use of NDVI_TOA or NDVI_RC can lead to the selection of the wrong K 443_670 and K 490_670 values (as shown in Tables 3-7) under some atmospheric conditions and can cause further errors in aerosol retrieval. Thus, to solve this problem, we developed a novel iterative AC method to calculate the real NDVI. Firstly, we used the Rayleigh scattering and gas correction NDVI_RC, which was calculated via Equation (16), as the initial NDVI value (recorded as NDV I init ); then, after combining the SCA and surface type, the initial K init 443_670 and K init 490_670 could be determined as: There were three problems to be solved in this step. The first was the determination of the solution τ for each aerosol model, the second was the determination of the NDVI, and the third was the determination of the best aerosol model. The steps are described in detail as follows: (1) For each aerosol model M and Θ, we found the solution τ(M, Θ) that made the cost function ε 1 (Equation (17)) equal to 0 (i.e., ε 1 = 0). Then, we obtained the R Cal SR (λ, M, Θ) of 443, 490, 670, and 865 nm: where R Cal SR (λ, M, Θ, τ) is the directional surface reflectance calculated using the AC method considering the non-Lambertian effect, as shown in step (c); the superscript i of NDV I i and K i 490_670 represents the number of iterations; and the NDV I init and K init 490_670 determined in step (d) are the initial values of the iteration (i.e., i = 0).
(3) We calculated the cost function ε 2 (Equation (19)) using R Cal SR (443, M, Θ), R Cal SR (670, M, Θ), and K i+1 430 670 M, NDV I i+1 , Θ , with some fitting error; the solution M was that which minimized ε 2 . Then, the corresponding multiangle mean value of τ(550 nm, Θ) was the retrieval solution τ(550 nm). Figure 6 shows the spatial distribution of the DPC AOD (550 nm) products retrieved from the VISRR on 3 April 2019 over eastern China; the main surface types of this region are grass, croplands, urban land, and water (lake and sea), as shown in Figure 1. Figure 6a is the DPC TOA true-color image (RGB) consisting of R Meas TOA (670), R Meas TOA (565), and R Meas TOA (490). It can be seen that most regions were cloud-free, except for some areas in the south and west, where an obvious strip of white cloud and some flocculent clouds existed; in addition, grey haze aerosols were distributed in most areas and blurred the surface. Figure 6b shows the DPC AOD (3.3 km) products: the blank areas represent those pixels identified as clouds, ice, and water; the blue regions with a low AOD (<0.2) indicate a clean atmosphere and were mainly distributed in the northern and eastern parts of the Shandong peninsula; the slight-pollution (0.25 < AOD < 0.5) regions were spread over southern Beijing-Tianjin-Hebei (BTH), western Shandong, northern Jiangsu, etc.; and moderate-and heavy-pollution regions with a high AOD (>0.5) could be found over most areas of Henan, central Anhui, and southern Jiangsu. Meanwhile, the DPC surface products' RGB consisting of R SR (670), R SR (565), and R SR (490) is shown in Figure 6c. It can be seen that the surface features became more distinguishable after removing the coupling effect of the atmosphere. In addition, the AOD products of the MYD04_3k_DT AOD (3 km), MYD04_L2_DB AOD (10 km), and MYD04_L2_DT AOD (10 km) results were chosen for comparison and are shown in Figure 6c-e, respectively. Overall, the spatial distribution between the DPC and MODIS AOD products was similar, but there were some differences. First, the AOD from the DPC and MYD04_L2_DB had high spatial coverage. However, the products from MODIS DT were missing in areas such as BTH and the Liaoning region, even if there was no cloud shown in the TOA reflectance RGB. The reason for this might be that the pixels were over-recognized as clouds, or that the DT algorithm cannot be implemented over bright land cover. Second, the spatial resolutions of the DPC and MYD04_3k_DT AOD were 3.3 km (individual pixel) and 3 km, respectively, while that of the MYD04_L2_DT and MYD04_L2_DB AOD was 10 km. In addition to the increased detail, the variations between different gradients of 3.3 or 3 km AOD products are also smoother than those of 10 km products. Third, the four AOD products were of different magnitudes on the whole; the DT AOD products had the highest value, while the DPC products had the lowest value. In this case, the DPC VISRR algorithm can simultaneously retrieve AOD products with the advantages of fine spatial resolution and high coverage.   Figure 7 shows another example of AOD products over western China and its surrounding areas on 30 May 2019, where the main surface types are barren land, grass, permanent snow, ice, and water (lake). From the TOA RGB (a), it can be seen that the permanent snow and ice cover was mainly located in the Himalayas, Kunlun, and the Tianshan Mountains. Thick, thin, and broken clouds were distributed in the Qinghai-Tibet Plateau and east of the Taklimakan Desert. Figure 7b presents the AOD products retrieved from the DPC, which show that the high aerosol loadings were distributed in the oases around the Taklimakan Desert and the Tarim River. Figure 7c-e display the AOD products of MODIS. It is evident that the DT algorithm is ineffective when applied to bright land cover compared to the VISRR and DB products. The DPC AOD products had a higher value over north Taklimakan and were more detailed than the DB products due to their fine spatial resolution.  Figure 7 shows another example of AOD products over western China and its surrounding areas on 30 May 2019, where the main surface types are barren land, grass, permanent snow, ice, and water (lake). From the TOA RGB (a), it can be seen that the permanent snow and ice cover was mainly located in the Himalayas, Kunlun, and the Tianshan Mountains. Thick, thin, and broken clouds were distributed in the Qinghai-Tibet Plateau and east of the Taklimakan Desert. Figure 7b presents the AOD products retrieved from the DPC, which show that the high aerosol loadings were distributed in the oases around the Taklimakan Desert and the Tarim River. Figure 7c-e display the AOD products of MODIS. It is evident that the DT algorithm is ineffective when applied to bright land cover compared to the VISRR and DB products. The DPC AOD products had a higher value over north Taklimakan and were more detailed than the DB products due to their fine spatial resolution.  Figure 8a shows the average DPC VISRR AOD products (March 2019 to March 2020) over China. There were four main high-value regions: (1) the northern and eastern plain of China, mainly including BTH, Shandong, Henan, Anhui, and Shanghai, where the aerosol emissions may be related to the fact that this region contains the largest population and the most industries; (2) the Sichuan Basin, mainly including eastern Sichuan and western Chongqing, where, in addition to industrial emissions, topographic and meteorological factors contribute to aerosol pollution [82,83]; (3) southern China, mainly including the south of Yunnan, Guangxi, and Hainan, where smoke aerosol pollution is transported in significant amounts from biomass burning in Southeast Asia [84,85]; (4) western China, mainly including the Taklimakan Desert and its surrounding areas-Taklimakan is the largest desert in China, and dust events often break out here, which are likely responsible for the aerosol pollution over this region [86]. Figure 8b-d present the AOD products from MYD04_3k_DT, MYD04_DB, and MYD04_DT and compare them with the MODIS products. In general, the spatial distributions of the AOD between the DPC and MODIS were similar. However, the AOD magnitude differed in that the MODIS AOD products were slightly higher than those of the DPC over northeastern China and the Qinghai-Tibet Plateau region. According to the validation of He et al. [47], MODIS 3 km and 10 km aerosol products are generally overestimated over China, and the error derives from two main sources. The first is the surface reflectance, and the second is the aerosol models in the LUT.

Validation Using Ground-Based Aerosol Products
In this study, twenty-two AERONET and SONET sites ( Figure 1 and Table 2) were chosen to validate the AOD (550 nm) products from the DPC and MODIS. Due to the lack of ground-based measurements with 550 nm AOD products, a second-order polynomial fitting method was chosen to calculate the 550 nm AOD using the 440, 500, and 675 nm values [87]. Furthermore, for the spatiotemporal matching between the ground and the satellites, the ground-based valid AOD products were averaged within one hour (30 min before or after the overpass of the satellite), and the satellite-based available AOD products were averaged within a sampling window (5 × 5 pixels) centered on the ground site [88]. Figure 9 shows the validation results for the DPC (a) and MODIS (MYD04_3k_DT (b), MYD04_L2_DB (c), and MYD04_L2_DB (d)) AOD products compared to the groundbased products from March 2019 to March 2020. y = a * x + b is a linear fitting equation, where a and b are the slope and intercept, respectively; N is the number of AOD pairs; and the good fraction (Gfrac) refers to the fraction of AOD pairs within the expected error EE = ±(0.05 + 20%). The evaluating indicators of the four aerosol products (a-d) were as follows: a and b were (0.96, 0.01), (1.03, 0.12), (0.98, 0.05), and (0.97, 0.1); the R was 0.88, 0.87, 0.82, and 0.87; the RMSE was 0.17, 0.21, 0.27, and 0.20; the N was 971, 1413, 2539, and 1997; and the Gfrac(s) was 62.7%, 51.7%, 60.6%, and 58.5%. The best evaluating indicators are presented in bold. The fact that the revisitation period of the DPC was longer than that of MODIS caused the total number of collected DPC AOD retrievals to be the lowest. At the same time, MYD04_L2_DB had the most retrieval results. The slopes of all products were close to 1; however, the DPC AOD products had the lowest intercept and RMSE and the highest R and Gfrac, demonstrating that the DPC AOD products had the highest accuracy over the study area. Specifically, compared with the AOD products at the 3 km level (i.e., Figure 9a

Validation Using Ground-Based Aerosol Products
In this study, twenty-two AERONET and SONET sites ( Figure 1 and Table 2) were chosen to validate the AOD (550 nm) products from the DPC and MODIS. Due to the lack of ground-based measurements with 550 nm AOD products, a second-order polynomial fitting method was chosen to calculate the 550 nm AOD using the 440, 500, and 675 nm values [87]. Furthermore, for the spatiotemporal matching between the ground and the satellites, the ground-based valid AOD products were averaged within one hour (30 min before or after the overpass of the satellite), and the satellite-based available AOD products were averaged within a sampling window (5 × 5 pixels) centered on the ground site [88]. The fact that the revisitation period of the DPC was longer than that of MODIS caused the total number of collected DPC AOD retrievals to be the lowest. At the same time, MYD04_L2_DB had the most retrieval results. The slopes of all products were close to 1; however, the DPC AOD products had the lowest intercept and RMSE and the highest R and Gfrac, demonstrating that the DPC AOD products had the highest accuracy over the study area. Specifically, compared with the AOD products at the 3 km level (i.e., Figure 9a versus Figure 9b), although they had a similar R value (0.88 vs. 0.87), the MYD04_3k_DT products had a low Gfrac (51.7%), a slope (1.03) greater than 1, and a high bias of intercept (0.12), revealing that these AOD products were generally overestimated. Compared with the AOD products at the 10 km level (i.e., Figure 9a versus Figure 9c,d), although they had similar Gfrac values (62.7%, 60.6%, and 58.5%), the MYD04_L2_DB products had a high RMSE (0.27) and a low R (0.82), and the MYD04_L2_DT products were similar to MYD04_3k_DT with a high bias of intercept (0.1). The slope smaller than 1 and the intercept greater than 0 implied that these AOD products were overestimated for small values but underestimated for large values. The validation results of three MODIS AOD products in this study were consistent with previous research [47,54]. high bias of intercept (0.12), revealing that these AOD products were generally overestimated. Compared with the AOD products at the 10 km level (i.e., Figure 9a versus Figure  9c,d), although they had similar Gfrac values (62.7 % , 60.6 %, and 58.5 % ), the MYD04_L2_DB products had a high RMSE (0.27) and a low R (0.82), and the MYD04_L2_DT products were similar to MYD04_3k_DT with a high bias of intercept (0.1). The slope smaller than 1 and the intercept greater than 0 implied that these AOD products were overestimated for small values but underestimated for large values. The validation results of three MODIS AOD products in this study were consistent with previous research [47,54]. To better understand and analyze the uncertainties of the DPC VISRR AOD products, Figure 10 shows the validation results over different land cover types: (a) urban, (b) cropland, (c) mixed, (d) forest, and (e) grassland. The evaluating indicators over the four surface types (a-d) are shown in Table 9, where the best evaluating indicators are presented in bold. There were 12 ground-based sites covering urban land; therefore, most of the collected DPC AOD retrievals were over this surface type. The slopes and intercepts of all the products were close to 1 and 0; the R was higher than 0.85, and the RMSE was lower than 0.22. Specifically, the accuracy of the DPC AOD over forests and croplands was highest among the surface types, with Gfrac values of 77.3% and 69.3%, respectively. The reasons for this can be summarized as follows: (1) compared to bright land cover, aerosol retrieval over vegetation (dark) surface types has advantages, as the surface signal accounts for a relatively low proportion of the VIS reflectance at TOA; (2) the VIS SRR with the NDVI and SCA constraints proposed in this study can well express the change in vegetation with time. The Gfrac (61.6%) of urban land was lower than that of forests and croplands, possibly because urban land cover is very complex, consisting of buildings, trees, water, roads, etc., and, with the rapid development of urbanization, urban land cover can change substantially over a short time [89]. Although the VISRR algorithm considers the discrepancies in the SRR over different surface types, it does not consider the change in the surface type over a short period of time, which can result in an estimation error for the surface that can be further transmitted to AOD retrieval. The mixed surface type was represented by three sites (Kashi, Jiaozuo, and Songshan): the Jiaozuo site was mainly a mix of urban land (56%) and cropland (46%); the Songshan site was mainly a To better understand and analyze the uncertainties of the DPC VISRR AOD products, Figure 10 shows the validation results over different land cover types: (a) urban, (b) cropland, (c) mixed, (d) forest, and (e) grassland. The evaluating indicators over the four surface types (a-d) are shown in Table 9, where the best evaluating indicators are presented in bold. There were 12 ground-based sites covering urban land; therefore, most of the collected DPC AOD retrievals were over this surface type. The slopes and intercepts of all the products were close to 1 and 0; the R was higher than 0.85, and the RMSE was lower than 0.22. Specifically, the accuracy of the DPC AOD over forests and croplands was highest among the surface types, with Gfrac values of 77.3% and 69.3%, respectively. The reasons for this can be summarized as follows: (1) compared to bright land cover, aerosol retrieval over vegetation (dark) surface types has advantages, as the surface signal accounts for a relatively low proportion of the VIS reflectance at TOA; (2) the VIS SRR with the NDVI and SCA constraints proposed in this study can well express the change in vegetation with time. The Gfrac (61.6%) of urban land was lower than that of forests and croplands, possibly because urban land cover is very complex, consisting of buildings, trees, water, roads, etc., and, with the rapid development of urbanization, urban land cover can change substantially over a short time [89]. Although the VISRR algorithm considers the discrepancies in the SRR over different surface types, it does not consider the change in the surface type over a short period of time, which can result in an estimation error for the surface that can be further transmitted to AOD retrieval. The mixed surface type was represented by three sites (Kashi, Jiaozuo, and Songshan): the Jiaozuo site was mainly a mix of urban land (56%) and cropland (46%); the Songshan site was mainly a mix of forests (36%), grassland (36%), and cropland (28%); and the Kashi site was mainly a mix of barren land (49%), cropland (40%), and grassland (20%). Over mixed-surface-type regions, there will be many mixed pixels. Although these pixels are identified as a certain surface type in advance, the selected prior SRRs may not represent the real relationship of the mixed pixels, which results in retrieval errors. Figure 10e shows the time series of the DPC AOD (blue dot) products against ground-based (red triangle) measurements over grassland (including the AOE_Baotou and Lhasa sites). The first 16 points are matchups from the AOE_Baotou site, and the rest are from the Lhasa site; the evaluating indicators N, RMSE, and Gfrac were 35, 0.1, and 62.9%, respectively. Although most of the valid matched AOD values were lower than 0.2 for AOD_Baotou and 0.1 for Lhasa, it can be seen that DPC AOD products were in better agreement with the ground-based measurements and showed the gradual process of AOD change.
Remote Sens. 2022, 14, x FOR PEER REVIEW 21 of 26 mix of forests (36%), grassland (36%), and cropland (28%); and the Kashi site was mainly a mix of barren land (49%), cropland (40%), and grassland (20%). Over mixed-surfacetype regions, there will be many mixed pixels. Although these pixels are identified as a certain surface type in advance, the selected prior SRRs may not represent the real relationship of the mixed pixels, which results in retrieval errors. Figure 10e shows the time series of the DPC AOD (blue dot) products against ground-based (red triangle) measurements over grassland (including the AOE_Baotou and Lhasa sites). The first 16 points are matchups from the AOE_Baotou site, and the rest are from the Lhasa site; the evaluating indicators N, RMSE, and Gfrac were 35, 0.1, and 62.9%, respectively. Although most of the valid matched AOD values were lower than 0.2 for AOD_Baotou and 0.1 for Lhasa, it can be seen that DPC AOD products were in better agreement with the ground-based measurements and showed the gradual process of AOD change.

Conclusions
In this study, an innovative algorithm called VISRR was proposed and applied to DPC multiangle intensity data to simultaneously retrieve AOD at an intrinsic spatial resolution (3.3 km) over dark and bright land cover. The VISRR algorithm considers surface anisotropy and neither requires a shortwave infrared band nor establishes SRDs. The VISRR's previous SRRs were built based on long-time-series MYD09 products and

Conclusions
In this study, an innovative algorithm called VISRR was proposed and applied to DPC multiangle intensity data to simultaneously retrieve AOD at an intrinsic spatial resolution (3.3 km) over dark and bright land cover. The VISRR algorithm considers surface anisotropy and neither requires a shortwave infrared band nor establishes SRDs. The VISRR's previous SRRs were built based on long-time-series MYD09 products and consideration of surface types, SCA, and NDVI. Meanwhile, to address the issue of the NDVI being susceptible to the atmosphere, an innovative iterative AC method was proposed to obtain realistic NDVI values. Then, by combining the fundamental aerosol models across China with a LUT method, the proposed algorithm was applied to thirteen months of DPC multiangle intensity data.
A long-time-series comparison and validation of the AOD products were performed, and the results showed that the VISRR algorithm is suitable for dark and bright land cover. The DPC and three MODIS AOD products had similar spatial distribution, but the DPC AOD products had both high spatial resolution and coverage. Validation using twenty-two ground-based sites over different surface types showed that the DPC AOD had the best performance (R of 0.88, RMSE of 0.17, and Gfrac of 62.7%) among the four products over the study area. Then, the uncertainties of the AOD products were comprehensively analyzed and discussed for future improvements. First, the DPC AOD products were validated over five typical surface types: forests had the best Gfrac (77.3%), while mixed land cover had the lowest Gfrac (53.5%). The possible source of error was the chance that the selected prior SRRs did not represent the real relationships of the mixed pixels.
In addition to the DPC/GF-5(01), VISRR can be applied to similar sensors, such as the advanced DPC/GF-5 (02), CM-1, DQ-1, DQ-2, and 3MI [49]. When combined with our previous fine-model aerosol algorithm SNOSPR, based on multiangle polarized data [30], the fine-model aerosol friction can be calculated, a key aerosol optical parameter retrieved from satellite remote sensing for distinguishing aerosol sources.