Retrieval and Validation of AOD from Himawari-8 Data over Bohai Rim Region, China

: The geostationary satellite Himawari-8, possessing the Advanced Himawari Imager (AHI), which features 16 spectral bands from the visible to infrared range, is suitable for aerosol observations. In this study, a new algorithm is introduced to retrieve aerosol optical depth (AOD) over land at a resolution of 2 km from the AHI level 1 data. Considering the anisotropic e ﬀ ects of complex surface structures over land, Moderate Resolution Imaging Spectroradiometer (MODIS) bidirectional reﬂectance distribution function (BRDF) model parameters product (MCD19A3) is used to calculate the surface reﬂectance for Himawari-8’s view angle and band. In addition, daily BRDF model parameters are calculated in areas with dense vegetation, considering the rapid variation of surface reﬂectance caused by vegetation growth. Moreover, aerosol models are constructed based on long duration Aerosol Robotic Network (AERONET) single scattering albedo (SSA) values to stand for aerosol types in the retrieval algorithm. The new algorithm is applied to AHI images over Bohai Rim region from 2018 and is evaluated using the newest AERONET version 3 AOD measurements and the latest MODIS collection 6.1 AOD products. The AOD retrievals from the new algorithm show good agreement with the AERONET AOD measurements, with a correlation coe ﬃ cient of 0.93 and root mean square error (RMSE) of 0.12. In addition, the new algorithm increases AOD retrievals and retrieval accuracy compared to the Japan Aerospace Exploration Agency (JAXA) aerosol products. The algorithm shows stable performance during di ﬀ erent seasons and times, which makes it possible for use in climate or diurnal aerosol variation studies.


Introduction
Atmospheric aerosols have profound impacts on the Earth's radiative balance and the climate through their direct radiative effects, indirect effects, and other cloud-mediated climatic effects [1,2]. On the other hand, as the primary pollutant, atmospheric aerosols have adverse impacts on people's health, since fine aerosol particles can be inhaled into human lungs, causing respiratory and cardiovascular diseases [3][4][5]. To study aerosol characteristics and their impacts to radiation, climate and human health, it is necessary to obtain in-depth observations and understand their physical and chemical characteristics.
Aerosol optical depth (AOD) is one of the important optical properties of aerosols, which describes the attenuation of sunlight by a total column of aerosols in the atmosphere [6]. Many satellite sensors are used them as inputs in the radiative transfer model to construct lookup tables for AOD retrieval. In order to solve the problem of surface reflectance changes due to the growth of vegetation, daily BRDF model parameters in areas with dense vegetation are calculated. In addition, average single scattering albedo (SSA) values are calculated under different air pollution conditions to constrain the aerosol model in the retrieval algorithm based on the long duration measurements from AERONET. The new AOD retrieval algorithm is applied to AHI level 1 data at a 2 km resolution during 2018. The retrieval algorithm is validated with AERONET AOD measurements and compared with the JAXA (Japan Aerospace Exploration Agency) and MODIS AOD products.

Study Region
The Bohai Rim region (about 518,000 km 2 ) is located in northern China (Figure 1), including the Hebei, Liaoning, and Shandong provinces, as well as the Beijing and Tianjin metropolitan areas. As one of the most populated regions in China, it has a population of 250 million. Due to the dense population and high industrialization, the region plays an important role in the development of China. In addition, Bohai Rim region has suffered from serious air pollution for many years.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 19 calculated based on MAIAC BRDF products, measurements of AHI TOA reflectance, and the Rossthick/Li-sparse (RTLS) model and are used them as inputs in the radiative transfer model to construct lookup tables for AOD retrieval. In order to solve the problem of surface reflectance changes due to the growth of vegetation, daily BRDF model parameters in areas with dense vegetation are calculated. In addition, average single scattering albedo (SSA) values are calculated under different air pollution conditions to constrain the aerosol model in the retrieval algorithm based on the long duration measurements from AERONET. The new AOD retrieval algorithm is applied to AHI level 1 data at a 2 km resolution during 2018. The retrieval algorithm is validated with AERONET AOD measurements and compared with the JAXA (Japan Aerospace Exploration Agency) and MODIS AOD products.

Study Region
The Bohai Rim region (about 518,000 km 2 ) is located in northern China (Figure 1), including the Hebei, Liaoning, and Shandong provinces, as well as the Beijing and Tianjin metropolitan areas. As one of the most populated regions in China, it has a population of 250 million. Due to the dense population and high industrialization, the region plays an important role in the development of China. In addition, Bohai Rim region has suffered from serious air pollution for many years.

AHI Data
Himawari-8, a geostationary satellite, was launched on the 7 October 2014 and is located at 140.7°E [34]. It carries an AHI including 16 spectral bands from visible to infrared with a spatial resolution of approximately 0.5 to 2 km ( Table 1). The temporal resolutions of the full-disk and target

AHI Data
Himawari-8, a geostationary satellite, was launched on the 7 October 2014 and is located at 140.7 • E [34]. It carries an AHI including 16 spectral bands from visible to infrared with a spatial resolution of approximately 0.5 to 2 km ( Table 1). The temporal resolutions of the full-disk and target area images are 10 and 2.5 min, respectively. Therefore, the Himawari-8/AHI can provide more detailed observations of atmospheric aerosols in Asia and Oceania, for instance, it is suitable for studying the diurnal variation of AOD [35]. The official JAXA P-Tree system provides the Himawari Standard Data (HSD format), level 1 gridded reflectance data (NetCDF4 format), and geophysical parameter data, which are readily accessible. In this study, daytime (10:00, 11:00, 12:00, 13:00, 14:00, and 15:00 in Beijing time) level 1 reflectance data at a spatial resolution of 2 km from 2018 are collected. The level 2 version 2.1 aerosol products at a 0.05 • (~5 km) spatial resolution were downloaded from the JAXA for comparison.

MODIS Data
The MODIS, loaded on the Terra and Aqua satellites, includes 36 spectral bands that extent from the visible to the infrared with spatial resolutions of 0.25, 0.5, and 1 km. The Terra and Aqua satellites cross the equator at about 10:30 a.m. and 1:30 p.m. (local time), respectively [36]. With a wide swath (~2330 km), the MODIS can observe the entire globe every 1 to 2 days.
Hereafter, "MXD" represents both MYD (Aqua) and MOD (Terra) products for convenience. The latest MODIS aerosol products (MXD04) were retrieved using the DT and DB algorithms with spatial resolutions of 10 and 3 km (3 km product only for DT). Some studies have shown that the AOD retrieved using the DB algorithm are generally more accurate than that by using the DT algorithm [10,11]. In addition, the DB algorithm can retrieve more efficient AOD over non-vegetated regions (i.e., urban and desert regions, etc.). So, MXD04 DB AOD products from January to December of 2018 were used to perform an intercomparison of AHI AOD retrievals in this study.
The new 8-day BRDF model parameter product MCD19A3 provides three weights (isotropic kernel, volumetric kernel, and geometric kernel) of the Ross-thick/Li-sparse (RTLS) BRDF model and can be used to describe the anisotropy of each pixel, also reporting the number of days since the last RTLS model update. The MCD19A3 product can be downloaded from the EarthData website (https://search.earthdata.nasa.gov/), so the data for 2016 and 2017 were used to build a new monthly surface anisotropic information database for AHI AOD retrieval in this study.

AERONET Data
The Aerosol Robotic Network (AERONET) project is a federation of ground-based remote sensing aerosol networks established by NASA and PHOTONS. It can provide accurate ground measurements of AOD every 15 min with an uncertainty of about 0.01 [37]. AERONET also provides other aerosol properties products, such as aerosol refractive index and single scattering albedo. The AERONET version 3 (V3) data were released in January 2018, and the AOD data are provided for three data quality levels: Level 1.0 (unscreened), level 1.5 (cloud-screened and quality controlled), and level 2.0 (quality-assured) [38]. In this study, the AERONET level 2.0 AOD data at three sites,

AOD Retrieval Methodology
Passive retrieval of AOD from satellites is usually based on the apparent reflectance ρ* at the top-of-atmosphere (TOA) as measured by a satellite sensor. Considering the surface BRDF effect, it can be expressed as follows [39,40]: where µ s and µ v are the cosines of the solar zenith angle and view zenith angle, respectively; φ is the relative azimuth angle; ρ 0 (µ s , µ v , φ) is the atmospheric path reflectance; τ is the optical depth of the atmosphere; T R+A (µ s ) and T R+A (µ v ) represent the downward and upward atmospheric transmittance, respectively; t d and S are the diffuse transmittance and the atmospheric spherical albedo, respectively; ρ t is the directional surface reflectance and ρ t , ρ t , and = ρ t represent the surface hemispherical-directional, directional-hemispherical, and hemispherical-hemispherical reflectance, respectively. Detailed descriptions of ρ t , ρ t and = ρ t can be found in Shi's work [41]. Equation (1) shows that the apparent reflectance at the TOA measured by the satellite sensor is attributed to the atmospheric path reflectance and the equivalent reflectance of the atmospheric and surface interactions. The atmospheric path reflectance includes the joint contribution of aerosols and molecules. Therefore, the aerosol model and LSR (calculated using the RTLS model based on prior knowledge of the BRDF) are crucial for satellite AOD retrieval over land. In addition, gas absorption and cloud detection also need to be carefully considered [42]. In this study, gas absorption is corrected using the midlatitude summer (winter) model defined in 6SV. Cloud detection is conducted based on the method proposed by Yang et al. with minor modification (4 pixels centered on the cloud pixel do not perform AOD retrieval) before further AOD retrieval [43].

Surface Reflectance Estimation
LSR is one of the key parameters for AOD retrieval over land. Previous studies have shown that an error of 0.01 in LSR leads to an error of 0.1 in AOD [44]. Methods using the lowest or second minimum value of MOD09 or TOA reflectance generally overestimate or underestimate surface reflectance due to anisotropic surface reflectance, as well as the possible significant changes of LSR due to the growth of vegetation [25,45].
Therefore, an a priori monthly BRDF database based on historical MCD19A3 data was built to calculate LSR for the AHI. Since the spectrum response function and wavelength range of the AHI band used is different from the MODIS sensor, spectral transformation was performed to minimize the effect of the different spectrum response between the two sensors.
In addition, considering that the LSR can significantly change due to the growth of vegetation in areas with dense vegetation, daily BRDF model parameters were calculated for AOD retrieval. The details of this LSR method are described below.

Method for Building Prior BRDF Database
In this study, the RTLS model is used to describe the directional reflectance of the land surface (i.e., the LSR). The model can be written as three terms: where ρ t (θ s , θ v , φ, λ) is the surface bidirectional reflectance; λ is the wavelength; f iso (λ), f vol (λ), and f geo (λ) represent the weight of isotropic scattering, radiative transfer-type volumetric scattering, Remote Sens. 2020, 12, 3425 6 of 19 and geometric-optical surface scattering, respectively; and K vol is the Ross-thick kernel, which can be calculated by Equation (3): where ξ is the phase angle and can be written as: K geo is the Li-sparse kernel, which can calculated by Equation (5): Notably, f iso (λ), f vol (λ), and f geo (λ) vary with time, but f iso (λ) change slowly, especially among the same period of years [32,46]. In addition, the BRDF shape (in the square bracket of Equation (2)) is related to the structure of the surface, while the reflectance magnitude (f iso (λ)) is determined by the microphysical properties of surface elements. Thus, in this study, we assumed that the BRDF shape and the reflectance magnitude remains unchanged among one month in areas with non-dense vegetation (i.e., urban areas and areas with sparse vegetation). The BRDF model parameters from the MODIS product, after a spectral transformation, can be used in these areas. In order to eliminate the possible outliers, the minimum standard deviation method was used to determine and f geo (λ) f iso (λ) of the BRDF based on 2 years of MCD19A3 data. The monthly average of f iso (λ) was used to present the monthly change in isotropic scattering. The minimum standard deviation method used in this study was operated as follows: (1) There were about 8 images (namely, 8 values for a pixel) for every month, and every set of 4 values was combined as a group after the invalid values were removed. (2) The combination with minimum standard deviation was selected and the mean value in the selected combination was used to determine the final value of each pixel.

Spectral Transformation
The spectral response functions and wavelengths used between the MODIS and AHI are different, which causes different LSR values to be found for the same band at the same pixel between the instruments. Therefore, spectral transformation was conducted. With the assumption that the LSR between the AHI and MODIS is linear, the LSR of the AHI can be expressed as follows: where ρ AHI λ and ρ MODIS λ are the actual LSR of a target at wavelength λ for the AHI and MODIS, respectively. The actual LSR can be calculated by Equation (10): where ρ i is the LSR in band i; Γ is the spectral response function; and S(λ) is the spectral curve. Spectral curves of 60 different features were collected from the ENVI spectral library. For both the AHI and MODIS, LSR values at the blue band were calculated. The LSRs for the AHI and MODIS were very close ( Figure 2) and the spectral convention model was then developed as follows: Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 19

Method of Daily BRDF Model Parameters Determination
Since LSR may significantly vary in regions with dense vegetation, using the monthly BRDF database will result in inaccurate AOD retrieval. Figure 3 shows the relationship between AOD and TOA reflectance in the blue band under different reflectance magnitudes. It indicates that a 0.02 error in the reflectance magnitude will lead to an error of ~0.10 in AOD retrieval when the AOD and reflectance magnitude are relatively low. The error will increase with the increase of the reflectance magnitude and decrease with the increase of AOD. Therefore, it is important to determine the BRDF parameters in a daily scale in areas with dense vegetation. The change of the LSR under the same observation conditions is mainly caused by the change of the microphysical properties of surface elements [32]. Thus, the BRDF shape from MODIS products can be used and the reflectance magnitude can be solved according to Equation (2) once the LSR is determined. Therefore, the key factor in determining the BRDF model parameters is to calculate the LSR in real time with acceptable accuracy. The BRDF model parameters will be used in the 6SV code to build the lookup tables (LUTs) for AOD retrieval. The DT algorithm demonstrates that there is an empirical relationship between the visible and near-infrared bands over regions with dense vegetation and that the LSR in blue band can be calculated based on AHI measurements. However, using the MODIS DT algorithm for AHI AOD retrieval will induce overestimation [27]. Recently, Wang developed optimal reflectance relationships between visible and SWIR bands with three NDVI classifications, illustrating improvement of the accuracy of AHI AOD retrieval [30]. Thus, we used the model built by Wang for areas with dense vegetation (NDVI > 0.33) to calculate the LSR in this study. Firstly, the LSR in the blue band was calculated based on the SWIR band reflectance for dense vegetation pixels. Secondly, the f iso λ was calculated according to Equation (2). Finally, the second lowest value of f iso λ during a day was selected as the daily value to remove the influence of possible aerosols/clouds on the calculation results.

Method of Daily BRDF Model Parameters Determination
Since LSR may significantly vary in regions with dense vegetation, using the monthly BRDF database will result in inaccurate AOD retrieval. Figure 3 shows the relationship between AOD and TOA reflectance in the blue band under different reflectance magnitudes. It indicates that a 0.02 error in the reflectance magnitude will lead to an error of~0.10 in AOD retrieval when the AOD and reflectance magnitude are relatively low. The error will increase with the increase of the reflectance magnitude and decrease with the increase of AOD. Therefore, it is important to determine the BRDF parameters in a daily scale in areas with dense vegetation. The change of the LSR under the same observation conditions is mainly caused by the change of the microphysical properties of surface elements [32]. Thus, the BRDF shape from MODIS products can be used and the reflectance magnitude can be solved according to Equation (2) once the LSR is determined. Therefore, the key factor in determining the BRDF model parameters is to calculate the LSR in real time with acceptable accuracy. The BRDF model parameters will be used in the 6SV code to build the lookup tables (LUTs) for AOD retrieval. The DT algorithm demonstrates that there is an empirical relationship between the visible and near-infrared bands over regions with dense vegetation and that the LSR in blue band can be calculated based on AHI measurements. However, using the MODIS DT algorithm for AHI AOD retrieval will induce overestimation [27]. Recently, Wang developed optimal reflectance relationships between visible and SWIR bands with three NDVI classifications, illustrating improvement of the accuracy of AHI AOD retrieval [30]. Thus, we used the model built by Wang for areas with dense vegetation (NDVI > 0.33) to calculate the LSR in this study. Firstly, the LSR in the blue band was calculated based on the SWIR band reflectance for dense vegetation pixels. Secondly, the f iso (λ) was calculated according to Equation (2). Finally, the second lowest value of f iso (λ) during a day was selected as the daily value to remove the influence of possible aerosols/clouds on the calculation results.

Aerosol Model
The aerosol model is another key parameter for AOD retrieval. An unsuitable assumption for the aerosol model will increase the uncertainty of AOD retrieval, potentially even leading to incorrect results. Continental aerosol model has been widely used to retrieve AOD in China [9,43,44]. However, the atmospheric conditions in Bohai Rim region are relatively complex, accompanied by heavy haze and dusty weather. The lookup table (LUT) built by the continental aerosol defined in the 6SV code can generate large errors in the retrieval of AOD at some periods in this area. In order to analyze the influence of the choice of aerosol model, the SSA, an important optical parameter of aerosols which reflects the scattering ability of aerosols, was set to different values to estimate the TOA reflectance with different AODs using the 6SV model. The relative errors of AOD retrieval with different values of the SSA are shown in Figure 4, in which the "true value" of the SSA was set as 0.93. It shows that the correct SSA information is important for AOD retrieval and that the retrieval errors increase with the increase of aerosol loading. In addition, when the SSA is underestimated, the error rises sharply. Thus, the observed SSA data from AERONET in Bohai Rim region during 2010-2018 were used to substitute the SSA (SSA of ~0.89 in the blue band) value of the continental aerosol model in the 6SV model to reduce the errors in the retrieval. In order to reflect the "dynamic" SSA (function of AOD), we separated all cases into two categories, either general (AOD < 0.8) or heavily polluted (AOD ≥ 0.8). The average SSA under the heavily polluted condition was calculated based on all the cases where the AOD was greater than or equal to 0.8. The average SSA under the general condition was calculated based on all the cases where the AOD was less than 0.8 but larger than 0.4 according to the recommendation from the AERONET team, since large SSA retrieval uncertainties may occur when the AOD is small. The two different average SSA values were used in the 6SV model for different atmospheric conditions (general or heavily polluted) for the construction of the LUTs. The results for the average SSA values in different bands are shown in Table 2.

Aerosol Model
The aerosol model is another key parameter for AOD retrieval. An unsuitable assumption for the aerosol model will increase the uncertainty of AOD retrieval, potentially even leading to incorrect results. Continental aerosol model has been widely used to retrieve AOD in China [9,43,44]. However, the atmospheric conditions in Bohai Rim region are relatively complex, accompanied by heavy haze and dusty weather. The lookup table (LUT) built by the continental aerosol defined in the 6SV code can generate large errors in the retrieval of AOD at some periods in this area. In order to analyze the influence of the choice of aerosol model, the SSA, an important optical parameter of aerosols which reflects the scattering ability of aerosols, was set to different values to estimate the TOA reflectance with different AODs using the 6SV model. The relative errors of AOD retrieval with different values of the SSA are shown in Figure 4, in which the "true value" of the SSA was set as 0.93. It shows that the correct SSA information is important for AOD retrieval and that the retrieval errors increase with the increase of aerosol loading. In addition, when the SSA is underestimated, the error rises sharply. Thus, the observed SSA data from AERONET in Bohai Rim region during 2010-2018 were used to substitute the SSA (SSA of~0.89 in the blue band) value of the continental aerosol model in the 6SV model to reduce the errors in the retrieval. In order to reflect the "dynamic" SSA (function of AOD), we separated all cases into two categories, either general (AOD < 0.8) or heavily polluted (AOD ≥ 0.8). The average SSA under the heavily polluted condition was calculated based on all the cases where the AOD was greater than or equal to 0.8. The average SSA under the general condition was calculated based on all the cases where the AOD was less than 0.8 but larger than 0.4 according to the recommendation from the AERONET team, since large SSA retrieval uncertainties may occur when the AOD is small. The two different average SSA values were used in the 6SV model for different atmospheric conditions (general or heavily polluted) for the construction of the LUTs. The results for the average SSA values in different bands are shown in Table 2.

Validation
For convenience, the AOD retrieved by the new method is denoted as AOD-new in this paper. The AERONET AOD measurements and the JAXA AOD product were used for validation and comparison purposes, and the MXD04 AOD products were used for intercomparison. The AHI AODs were averaged within a 3 × 3 pixel grid centered around each AERONET site. The AERONET AODs with a time difference between the AHI AODs less than 15 min were averaged and interpolated to 550 nm to match the retrieved AOD wavelength [10]. Since the JAXA AODs were retrieved at the wavelength of 500 nm, the AERONET data during 2015-2017 were used to build a spectral transformation model. The result ( Figure 5) shows that the simple linear transformation model performs well.

Validation
For convenience, the AOD retrieved by the new method is denoted as AOD-new in this paper. The AERONET AOD measurements and the JAXA AOD product were used for validation and comparison purposes, and the MXD04 AOD products were used for intercomparison. The AHI AODs were averaged within a 3 × 3 pixel grid centered around each AERONET site. The AERONET AODs with a time difference between the AHI AODs less than 15 min were averaged and interpolated to 550 nm to match the retrieved AOD wavelength [10]. Since the JAXA AODs were retrieved at the wavelength of 500 nm, the AERONET data during 2015-2017 were used to build a spectral transformation model. The result ( Figure 5) shows that the simple linear transformation model performs well.

Overall Validation
Firstly, the AOD values retrieved by the new algorithm were compared against the AERONET AOD measurements using all 1956 matched samples at three sites during 2018. The correlations between the satellite-retrieved AODs (AOD-new and JAXA AOD) and AERONET AODs are shown in Figure 6a The bias between the satellite-retrieved AODs (AOD-new and JAXA AOD) versus AERONET AODs, varying with the AOD, is shown in Figure 7, in which all matchup samples have been divided into 10 bins with 200 samples in each bin (156 samples in the 10th bin). As shown in the Figure 7a, the bias of the AOD-new decreases with increase of the AERONET AOD when the AERONET AOD is smaller than 0.5. Figure 7b shows the bias between the JAXA AOD versus the AERONET AOD. Large biases and uncertainties are found at all aerosol loadings situation. Both the AOD-new and JAXA AODs show overestimation when the atmosphere is clean, but the bias of AOD-new is much smaller than that of JAXA AOD. Figure 7c,d is the same as Figure 7a,b, but varying with satelliteretrieved AODs. Figure 7c shows that the bias of AOD-new does not significantly vary with the increase of the AHI AOD and that there is slight overestimation over the area with high retrieved value. Compared with the bias of AOD-new, the bias of JAXA AOD, shown in Figure 7d, varies more dramatically, almost showing a linear increase trend with the increase of the JAXA AOD. The bias between the satellite-retrieved AODs (AOD-new and JAXA AOD) versus AERONET AODs, varying with the AOD, is shown in Figure 7, in which all matchup samples have been divided into 10 bins with 200 samples in each bin (156 samples in the 10th bin). As shown in the Figure 7a, the bias of the AOD-new decreases with increase of the AERONET AOD when the AERONET AOD is smaller than 0.5. Figure 7b shows the bias between the JAXA AOD versus the AERONET AOD. Large biases and uncertainties are found at all aerosol loadings situation. Both the AOD-new and JAXA AODs show overestimation when the atmosphere is clean, but the bias of AOD-new is much smaller than that of JAXA AOD. Figure 7c,d is the same as Figure 7a,b, but varying with satellite-retrieved AODs. Figure 7c shows that the bias of AOD-new does not significantly vary with the increase of the AHI AOD and that there is slight overestimation over the area with high retrieved value. JAXA AODs show overestimation when the atmosphere is clean, but the bias of AOD-new is much smaller than that of JAXA AOD. Figure 7c,d is the same as Figure 7a,b, but varying with satelliteretrieved AODs. Figure 7c shows that the bias of AOD-new does not significantly vary with the increase of the AHI AOD and that there is slight overestimation over the area with high retrieved value. Compared with the bias of AOD-new, the bias of JAXA AOD, shown in Figure 7d, varies more dramatically, almost showing a linear increase trend with the increase of the JAXA AOD.

Validation at Different Sites
The performances of the AHI AOD products (AOD-new and JAXA AOD) were examined at three sites (Table 3). Beijing and Beijing-CAMS are typical urban sites, while Xianghe is a countryside site near croplands. The comparison results show high correlation between the AOD-new and AERONET AOD data at all the three sites, with R values of 0.91, 0.92, and 0.94 and RMSE values of 0.12, 0.12, and 0.12, respectively. However, the best agreement of the JAXA AOD to the AERONET AOD was found at the Beijing-CAMS site, with the R of 0.73 and the RMSE of 0.32. The new retrieval algorithm achieved more than 50% decrease for the RMSE, MAE, and MRE compared to the JAXA method. In addition, the JAXA AOD exhibits the worst performance at the Xianghe site, with the R of 0.66 and RMSE of 0.42. In summary, these results demonstrate that the new method performs well over both bright and dark surface sites when utilizing the Himawari-8 data.

Validation in Different seasons
The performances of the AOD-new and JAXA AOD data in different seasons are shown in Figure 8. The correlation coefficients between the AOD-new (JAXA) and AERONET AOD data in spring, summer, autumn, and winter were 0.91 (0.79), 0.97 (0.90), 0.96 (0.78), and 0.86 (0.38), respectively. The new algorithm performs much better than the JAXA algorithm during all seasons. The RMSE, MAE, and MRE values for the new algorithm are smaller and the RMB values of AOD-new are closer to 1. It is noticed that there are some samples overestimation or underestimation largely in spring when aerosol loading is heavy. This indicates that the aerosol type changed significantly in this season. The JAXA AOD showed the best performance during summer, with the RMSE of 0.16 and MAE of 0.12, respectively. However, less agreement was found during autumn and winter. The RMSE, MAE, and RMB values in winter were 0.42, 0.30 and 2.08, respectively. These results demonstrate that the AHI AOD retrieved using the new algorithm feature considerable accuracy during each season, while the JAXA AOD performs poorly in autumn and winter.

. Validation in Different seasons
The performances of the AOD-new and JAXA AOD data in different seasons are shown in Figure  8. The correlation coefficients between the AOD-new (JAXA) and AERONET AOD data in spring, summer, autumn, and winter were 0.91 (0.79), 0.97 (0.90), 0.96 (0.78), and 0.86 (0.38), respectively. The new algorithm performs much better than the JAXA algorithm during all seasons. The RMSE, MAE, and MRE values for the new algorithm are smaller and the RMB values of AOD-new are closer to 1. It is noticed that there are some samples overestimation or underestimation largely in spring when aerosol loading is heavy. This indicates that the aerosol type changed significantly in this season. The JAXA AOD showed the best performance during summer, with the RMSE of 0.16 and MAE of 0.12, respectively. However, less agreement was found during autumn and winter. The RMSE, MAE, and RMB values in winter were 0.42, 0.30 and 2.08, respectively. These results demonstrate that the AHI AOD retrieved using the new algorithm feature considerable accuracy during each season, while the JAXA AOD performs poorly in autumn and winter.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 19   Figure 9 shows the diurnal variation of AOD and the evaluation metrics, which were calculated using all matched samples at three sites during 2018. Figure 9a shows that AOD retrieved by the new method are close to AERONET measurements, with little overestimation or underestimation at different hours in a day. The aerosol loadings decrease from 10:00 (Beijing time) to 12:00, and there is an obviously increase from 14:00 to 15:00. Our AOD retrievals have the same time trend with the AERONET measurements, while the JAXA AOD product presents an opposite trend. In addition, there is an obvious bias between the JAXA and AERONET AOD data. According Figure 9b, the retrieval accuracy of our method is stable for different times. There are high correlations between the AOD-new and AERONET measurements, with R values varying from 0.91 (10:00) to 0.95 (14:00). The MAEs and RMSEs show a decreasing trend from 10:00 to 15:00. The evaluation metrics of the JAXA  Figure 9 shows the diurnal variation of AOD and the evaluation metrics, which were calculated using all matched samples at three sites during 2018. Figure 9a shows that AOD retrieved by the new method are close to AERONET measurements, with little overestimation or underestimation at different hours in a day. The aerosol loadings decrease from 10:00 (Beijing time) to 12:00, and there is an obviously increase from 14:00 to 15:00. Our AOD retrievals have the same time trend with the AERONET measurements, while the JAXA AOD product presents an opposite trend. In addition, there is an obvious bias between the JAXA and AERONET AOD data. According Figure 9b, the retrieval accuracy of our method is stable for different times. There are high correlations between the AOD-new and AERONET measurements, with R values varying from 0.91 (10:00) to 0.95 (14:00). The MAEs and RMSEs show a decreasing trend from 10:00 to 15:00. The evaluation metrics of the JAXA validation results show that the JAXA method performs worse during the afternoon. The difference in accuracy for the JAXA AOD product at different times is mainly caused by changes in the geometric angle, which can affect the direction reflectance. As shown in Figure 10a, the scattering angle decreases between 10:00 and 15:00. In addition, the evaluation metrics shown in Figure 10b demonstrate that neither method performs very well when the scattering angle is high; however, our results still show a clear advantage. These results suggest that the new method can better represent the LSR values at different angles. Therefore, the method can catch the diurnal variation of AOD well when compared to the JAXA method. The spatial distribution of the retrieved hourly AOD is shown in Figure A1 in Appendix A.

Validation at Different Times
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 19 validation results show that the JAXA method performs worse during the afternoon. The difference in accuracy for the JAXA AOD product at different times is mainly caused by changes in the geometric angle, which can affect the direction reflectance. As shown in Figure 10a, the scattering angle decreases between 10:00 and 15:00. In addition, the evaluation metrics shown in Figure 10b demonstrate that neither method performs very well when the scattering angle is high; however, our results still show a clear advantage. These results suggest that the new method can better represent the LSR values at different angles. Therefore, the method can catch the diurnal variation of AOD well when compared to the JAXA method. The spatial distribution of the retrieved hourly AOD is shown in Figure A1 in Appendix A.     Figure 11 shows the spatial distribution of AOD from different retrieval algorithms on 1 November 2018. The first column shows the AHI true color images. The second to last column are images of AOD from AOD-new, JAXA AOD and the MODIS DB AOD (MXD04) data, respectively. The spatial resolutions of the AOD-new, JAXA AOD, and MXD04 AOD data are 2, 5, and 10 km, respectively. Compared with the JAXA AOD, the available retrievals and spatial distribution of the AOD-new are closer to the MXD04 DB AOD. In contrast, the JAXA AOD shows a large amount of missing data due to the retrieval algorithm that is used. According to Yoshida's work [25], land surface reflectance is estimated mainly by using the second lowest reflectance at 470 nm within a month, which may introduce noise that may cause retrieval failure. There are significant differences in the distribution between the JAXA AOD product and the other two data at 13:00. It is shown that the values of the JAXA AOD are much higher than both the AOD-new and MXD04 AOD. In addition,  Figure 11 shows the spatial distribution of AOD from different retrieval algorithms on 1 November 2018. The first column shows the AHI true color images. The second to last column are images of AOD from AOD-new, JAXA AOD and the MODIS DB AOD (MXD04) data, respectively. The spatial resolutions of the AOD-new, JAXA AOD, and MXD04 AOD data are 2, 5, and 10 km, respectively. Compared with the JAXA AOD, the available retrievals and spatial distribution of the AOD-new are closer to the MXD04 DB AOD. In contrast, the JAXA AOD shows a large amount of missing data due to the retrieval algorithm that is used. According to Yoshida's work [25], land surface reflectance is estimated mainly by using the second lowest reflectance at 470 nm within a month, which may introduce noise that may cause retrieval failure. There are significant differences in the distribution between the JAXA AOD product and the other two data at 13:00. It is shown that the values of the JAXA AOD are much higher than both the AOD-new and MXD04 AOD. In addition, with a higher spatial resolution, the AOD-new can reflect more detailed information for AOD distributions. with a higher spatial resolution, the AOD-new can reflect more detailed information for AOD distributions.     Figure 12b, the spatial distributions of the JAXA AODs and MXD AODs show many more differences, with the R of 0.60 and the RMSE of 0.53. Figure 12a,b show scatter plots of the AOD-new versus MXD04 AOD and JAXA AOD versus MXD04 AOD, respectively. The results from 548,283 matched samples show good agreement between the AOD-new and MXD04 AOD, with the R of 0.91. As seen in the Figure 12a, most of the samples are concentrated around the 1:1 line. The RMSE, MAE, and RMB values between the AODnew and MXD04 AOD are 0.18, 0.12, and 0.96, respectively. From Figure 12b, the spatial distributions of the JAXA AODs and MXD AODs show many more differences, with the R of 0.60 and the RMSE of 0.53.

Discussions
The LSR and aerosol model are key parameters for AOD retrieval. In this study, a new method considering the effect of the surface BRDF and local aerosol type is presented. Compared to the methods under the Lambertian surface approximation (such as the DT method, minimum albedo/reflectance method, and some modified DT methods), the new method is more realistic because the anisotropic reflectance of the land surface is considered. Compared with studies that simultaneously retrieve the BRDF (or direction reflectance) and AOD, we calculated the BRDF model

Discussion
The LSR and aerosol model are key parameters for AOD retrieval. In this study, a new method considering the effect of the surface BRDF and local aerosol type is presented. Compared to the methods under the Lambertian surface approximation (such as the DT method, minimum albedo/reflectance method, and some modified DT methods), the new method is more realistic because the anisotropic reflectance of the land surface is considered. Compared with studies that simultaneously retrieve the BRDF (or direction reflectance) and AOD, we calculated the BRDF model parameters independently before determining the AOD to reduce the computation burden. Instead of using the monthly BRDF database in all areas [33], daily BRDF model parameters were calculated according to short-wave infrared (SWIR) data and the RTLS model in areas with dense vegetation to solve the problem of surface reflectance changes due to the growth of vegetation. For aerosol model determination, we determined SSA values with the AOD according to the AERONET SSA at different AOD ranges in the study area. This is more suitable than selecting an aerosol model by changing the AOD and aerosol models to find the minimum errors between measurements and forward simulations of AHI bands. This is because the candidate aerosol models used as the final aerosol model may be unable to represent the real local aerosol type, which will lead to retrieval errors in practice.
The validation of AOD retrievals against AERONET AOD measurements, shown in Section 4, shows that the new method performs better than the JAXA algorithm (with R, RMSE, and MAE values of 0.93, 0.12, and 0.08, respectively, while the related metrics of the JAXA algorithm are 0.69, 0.37 and 0.25, respectively). In addition, it performs robustly with different sites, seasons, and times. However, increasing errors are shown when the scattering angle exceeds 160 • . This is because the LSR is higher when scattering is backward, and according to the study of Sun et al. [44], errors of AOD retrieval caused by LSR errors will increase when the LSR increases. There are few ground sites in the study area and most of the regions cannot be validated. Thus, MODIS aerosol products have been used for intercomparison, and the results show good agreement between the AOD-new and MXD04 AOD.

Conclusions
A new AOD retrieval algorithm, using prior knowledge of the BRDF, has been proposed to improve Himawari-8 AOD retrieval in this study. To reduce the effect of surface anisotropy on AOD retrieval, a BRDF parameter database was built using a minimum standard deviation synthesis method based on historical MCD19A3 data. A daily BRDF model parameter determination method was developed for areas with dense vegetation. Meanwhile, the continental aerosol model defined in the 6SV code was modified by the AERONET SSA parameters to obtain more accurate LUTs. Based on the new method, AOD were retrieved from Himawari-8 images with a spatial resolution of 2 km.
The validation results suggest that the AOD values retrieved by new method have good agreement with the AERONET AOD measurements, with an overall correlation coefficient of 0.93. The accuracy at each ground site is high, no matter whether the site is at bright or dark surface. The validation results also show that the new method performs well during all seasons and that the accuracy is not sensitive to time during the day. Through intercomparison with two other aerosol products, we have found that the retrieved AOD data provide a wider spatial coverage and can reflect more detailed information than the JAXA AOD product. The new AOD retrieval has high correlation with the MODIS AOD product. These results demonstrate that MODIS BRDF products can be used for AHI AOD retrieval, even when there are differences in the observation angles between two sensors. In addition, accurate AOD retrieval can improve PM 2.5 (particulate matter with an aerodynamic diameter < 2.5 µm) estimation, air quality forecasting models, and, additionally, diurnal aerosol variation studies. However, there are still some challenges in the retrieval, especially when the scattering angle approaches 180 degrees.