Preliminary Investigation of a New AHI Aerosol Optical Depth ( AOD ) Retrieval Algorithm and Evaluation with Multiple Source AOD Measurements in China

The Himawari-8 geostationary weather satellite, which is an Earth observing satellite launched in October 2014, has been applied in climate, environment, and air quality studies. Using hourly observation data from the Advanced Himawari Imager (AHI) on board Himawari-8, a new dark target algorithm was proposed to retrieve the aerosol optical depth (AOD) at 1 km and 5 km resolutions over mainland China. Because of the short satellite operation time and lack of AErosol RObotic NETwork (AERONET) sites across China, we cannot derive robust and representative surface reflectance relationships for the visible to near-infrared channels by atmospheric correction. Therefore, we inherited the empirical reflectance relationship from the Moderate Resolution Imaging Spectroradiometer (MODIS) and we used the AHI and MODIS spectral response functions to make the relationship more suitable for AHI. Ultimately, our AOD products can better reflect the regional characteristics with the AHI sensor. Seasonal averages showed that our product is more similar to MODIS Collection 6 (C6) Dark Target (DT) AOD than the Japan Aerospace Exploration Agency (JAXA) AHI AOD, but the difference is largest in winter. In addition, we evaluated several satellite retrieval products (our AHI AOD, JAXA AHI AOD, the National Oceanic and Atmospheric Administration (NOAA) VIIRS AOD, MODIS DT AOD, and MODIS DB AOD) against AERONET AOD from July 2016 to June 2017. The results showed that our AHI measurements demonstrate good agreement with, but exhibit a little overestimation, as compared to ground-based AERONET measurements with a correlation coefficient of 0.83 and an root-mean-square error (RMSE) of 0.112. The hourly validation also showed stable statistical results. A time series comparison with ground-based observations from two AERONET sites (Beijing-CAMS and XiangHe) showed that our AHI AOD products have trends as those in MODIS DB AOD, but that the bias in Beijing-CAMS is positive and higher than that in XiangHe. This error and the slight overestimation may be caused by the single continental aerosol model assumption and not considering the Normalized Difference Vegetation Index (NDVI).


Introduction
Aerosols are liquid or solid particles that are suspended in the atmosphere and are the most variable and uncertain component of the air.Generic aerosol types (urban and industrial regional pollution, smoke from vegetation fires, desert dust blown into the atmosphere by wind, and oceanic aerosol) are associated with different sources and emission mechanisms that are expected to exhibit optically significant differences [1,2].Aerosols are intricately linked to the climate system and the hydrologic cycle [2], and force the global energy budget [3], perturb the hydrological cycle [4], and are detrimental to human health in large concentrations [5].In addition, aerosol particles play a role in many chemical reactions in the atmosphere.If a large number of particles (especially toxic particles) are suspended in the atmosphere, it will have an impact on air quality.However, the scientific understanding of roles of aerosols in climate change is poor, and large uncertainties exist.Satellite observation is an advanced technology that makes important contributions to research in these areas [6].Furthermore, due to the heterogeneous distribution of sources, short lifetime, and episodic features of emission events, aerosols exhibit high spatiotemporal variability.Therefore, continuous global and regional characterization can only be realized through satellite remote sensing.
Remote sensing has been increasingly used to obtain aerosol properties, such as aerosol optical depth (AOD) that was researched in this work over the past decades [7].Although there is some uncertainty (sensor differences, calibration/characterization differences, retrieval algorithm differences, pixel selection, cloud, and other masking, etc.) in satellite retrievals, it still continues to develop [8][9][10][11][12].Since the first operational aerosol products were generated using data from the radiometer aboard the Tiros-N satellite that was launched on 19 October 1978, several instruments have observed aerosols from space [13].The Advanced Very High Resolution Radiometer (AVHRR) series of instruments were initially used only to provide AOD over oceans, but algorithms have been developed to provide AOD over land [14][15][16] for more than 35 years.The dual view Along-Track Scanning Radiometers (ATSR-2) and AATSR provided AOD over land and oceans from 1995 to 2012 [17,18].These time series data were extended with the Sea and Land Surface Temperature Radiometer (SLSTR), which was launched in 2016.The Moderate Resolution Imaging Spectroradiometer (MODIS) on board the Earth Observing System (EOS) satellite sensors provided the first operational global satellite dataset over land [6,19] using two instruments that were launched on the Terra satellite in December 1999 and on the Aqua satellite in May 2002.Of course, the uncertainty of the MODIS Collections (Collection 6 (C6), Collection 5.1 (C5.1)) aerosol datasets and the differences between the MODIS C6 Dark Target (DT) and C6 Deep Blue (DB) retrievals have been mentioned [8,11].Combining the AOD from the ATSR and MODIS/Terra instruments provides an AOD time series of more than 20 years [20].More recently, the retrievals were made with the Visible Infrared Imaging Radiometer Suite (VIIRS) on board the Suomi National Polar orbiting Partnership (S-NPP) satellite [21,22], which extended the MODIS time series [23].
In addition to polar orbiting satellites, aerosol retrievals have also been made from imaging instruments on geostationary platforms.The GOES Aerosol/Smoke Product (GASP) is a retrieval of AOD using visible imagery and is currently being used to monitor particulate pollution at a high temporal resolution.Wang et al. [24] used GOES-8 imagery to retrieve dust aerosol optical thickness over the Atlantic Ocean during the Puerto Rico Dust Experiment (PRIDE).In addition, Prados et al. [25] evaluated the GASP AOD from the GOES-12 imager over North America at various temporal and spatial scales based on comparisons with AOD from the AErosol RObotic NETwork (AERONET) and MODIS.Govaert et al. [26] estimated the surface albedo increase during the Sahel drought in the 1980s from Meteosat observations.Frequent observations from geostationary satellites are important in observing the diurnal evolution of aerosols and their effect on air quality, but polar orbiting satellites have insufficient temporal resolution to fully quantify these aerosol effects [27], and GOES cannot view the Asia-Pacific region [28].Conversely, the Advanced Himawari Imager (AHI) on board the Japanese Himawari-8 geostationary weather satellite provides an opportunity to observe the diurnal changes in aerosol properties over the Asia-Pacific region.In this paper, we present a method and demonstrate AOD retrieval using AHI data.
AOD retrieval is influenced by different factors, including calibration accuracy, the presence of clouds, the reflectance of the underlying surface, and aerosol properties [29].In this work, a spatial variability method inherited from MODIS was used to mask clouds.In addition, a universal aerosol mode, the continental aerosol mode, was set as the background aerosol conditions.However, it is extremely difficult to infer the part of the signal that is related to aerosols from the total signal that was detected by a satellite instrument, which is composed of contributions from both the atmosphere and the underlying surfaces.Especially affected by the adjacency effect that is caused by the scattering of reflected radiance from an inhomogeneous surface, it can introduce a blurring effect on remotely sensed data [30].Our algorithm that is described in this paper is based on the DT algorithm, which assumes that the surface reflectances in the visible and near-infrared bands exhibit a stable relationship [31].The empirical relationship for the MODIS data was calculated using over 15,000 pairs of MODIS/AERONET data at over 200 global sites [32].However, because of the short AHI operation time and the lack of AERONET sites over China, we cannot derive the correct relationship directly using AHI/AERONET collocations.To overcome this problem, we inherited the MODIS relationship using a channel transformation method.Based on the AHI reflectance in three channels (visible channels: 0.46 µm and 0.64 µm and near-infrared channel: 2.30 µm), we used a DT algorithm to retrieve AOD over land in China.Then, this work evaluates the new AHI AOD from July 2016 to June 2017 at different temporal and spatial scales based on comparisons with AOD values from different sources.
In this paper, we made a continental aerosol model assumption, used a cloud mask method based on spatial variability, and inherited the surface reflectance relationship from MODIS to develop a preliminary AOD retrieval algorithm for AHI.Section 2 introduces relevant satellite data and the AOD product, and Section 3 introduces the main methodology.The spatial and temporal comparison of AHI AOD with different sources and the validation against ground-based AERONET data are described in Section 4. Finally, we discuss the deficiencies of our AOD product and conclude the paper in Section 5.

Himawari-8/the Advanced Himawari Imager (AHI)
Himawari-8 is a Japanese weather satellite that is operated by the Japan Meteorological Agency (JMA) and was launched on 17 October 2014.The AHI on board Himawari-8, which was designed and built by Exelis Geospatial Systems (now Harris Space & Intelligence Systems), is the primary instrument.The AHI, at a spatial resolution of 0.5-2 km and a temporal resolution of 10-2.5 min, is a 16-channel (Table 1) multispectral imager that is intended to capture the visible light and infrared images over the Asia-Pacific region for 15 years.The AHI has similar spectral and spatial characteristics as the Advanced Baseline Imager (ABI) that is planned for use in the American GOES-R satellites.For this work, we use the AHI visible bands (0.46 µm, 0.64 µm) and near-infrared band (2.30 µm) for aerosol retrieval (1 km and 5 km resolution AOD product) and use Himawari-8 level three hourly AOD data from the Japan Aerospace Exploration Agency (JAXA) for comparison with our retrieval results.The AHI AOD products are 5 km grid data.For quality assurance, the AODs have four confidence levels: "very good", "good", "marginal", and "no confidence" (or "no retrieval").We obtained free data from the following website: (http://www.eorc.jaxa.jp/ptree/index.html).

The Moderate Resolution Imaging Spectroradiometer (MODIS) and the Visible Infrared Imaging Radiometer Suite (VIIRS)
The MODIS instrument [33] flies on EOS Terra (on a descending orbit (southward) over the equator at approximately 10:30 local sun time) and Aqua (on an ascending orbit (northward) over the equator at approximately 13:30 local sun time) polar orbiting satellites.The MODIS instruments measure spectral radiance in 36 channels (between 0.412 to 14.2 µm) in resolutions of between 250 m and 1 km (at nadir).From a vantage of approximately 700 km above the Earth and a with ±55 • view scan, each MODIS views a swath approximately 2330 km, and thus, nearly the entire globe is observed at least once daily, and the return period is 16 days.
VIIRS was launched aboard the S-NPP satellite in October 2011 and it is sensitive to 22 wavelength bands ranging in central wavelength from 0.412 µm to 12.01 µm [21].The instrument has a field of view of 112.56 • in the across-track direction and a nominal altitude of 824 km, corresponding to a swath width of approximately 3060 km with no data gaps near the equator.S-NPP was launched into sun-synchronous orbit with an equatorial crossing time of 13:30 (local time).
For this work, we made extensive use of MODIS/AQUA L2 (https://ladsweb.modaps.eosdis.nasa.gov/search/) and NOAA VIIRS-EDR (https://www.class.ncdc.noaa.gov/saa/products/welcome)products.All MODIS L2 (MOD04 L2 for Terra and MYD04 L2 for Aqua, which is known commonly as M*D04) aerosol products from two different algorithms are provided over land, known as DB and DT [19,34,35].The environmental data record (EDR) is the official level 2 product of VIIRS, and the AOD product also uses DT algorithm.We choose data with different quality flags (MODIS DB, MODIS DT, and VIIRS AOD, all quality data, and very good data) in this work, and more detail is listed in Table 2 [6,23].

The AErosol RObotic NETwork
The sun photometers of the AERONET program [36][37][38] constitute a network of ground-based remote sensing aerosol networks that were established by the National Aeronautics and Space Administration (NASA) and PHOTONS (PHOtométrie pour le Traitement Opérationnel de Normalisation Satellitaire).The AERONET collaboration can provide globally distributed observations of spectral AOD computed at three data quality levels: level 1.0 (unscreened), level 1.5 (cloud screened), and level 2.0 (cloud screened and quality assured) [32].Level 2.0 AERONET data were chosen for this study.The utilized AERONET optical thickness (τ) values are calculated by taking the mean of the daytime AERONET retrievals half an hour before and after satellite transit time (at least four observations).If for any reason there are no available AERONET τ values, the corresponding VIIRS/MODIS/Himawari-8 overpass is not included in this study.
Because the AHI operates at a wavelength of 0.51 µm, AERONET-observed values are interpolated to 0.51 µm using the Ångström exponent (AE) relationship after solving for α (Equation ( 1)) [39,40].The Ångström exponent, which is used to describe the dependency of aerosol optical thickness or aerosol extinction coefficient on wavelength, is inversely related to the average aerosol particle size.
where α is the Ångström exponent.τ λ 1 and τ λ 2 are the measurements of optical thickness at two different wavelengths and λ 1 and λ 2 , respectively.To obtain τ λ , Equation ( 1) is rearranged such that: where τ λ is the optical thickness at wavelength λ and τ λ 0 is the optical thickness at the reference wavelength λ 0 .
In this paper, all satellite data should match the in situ observations from July 2016 to June 2017 in the areas around China (including Taiwan and Hong Kong).Within a 25-km longitude and latitude scope, more than 25% of the total number pixels are used to calculate the mean.Using the value of the AERONET site within half an hour before and after the satellite transit and for values with a corresponding number of ground-based site data points that is no less than 4, we calculate the mean as a transit time site value.
Hereafter, for convenience, the Aqua/MODIS DT AOD is denoted as AOD_DT, the Aqua/MODIS DB AOD is denoted as AOD_DB, the NOAA S-NPP/VIIRS AOD is denoted as AOD_V, the JAXA H8/AHI AOD is denoted as AOD_JAXA, our new H8/AHI AOD is denoted as AOD_NEW, and the AERONET AOD at 550 nm is denoted as AOD_AERO.

Dark Target (DT) Algorithm
Kaufman and Tanré [41] first proposed a method for retrieving aerosol data over dense dark vegetation (known as the DT algorithm) using a multiband algorithm from the MODIS satellite images.Because of liquid water absorption and chlorophyll reflectance, there is an empirical relationship between the visible and the near-infrared channels [42].This algorithm assumes the land surface to be Lambertian and utilizes the surface reflectance relationship to decouple the atmospheric signal from the combined surface/atmosphere signal.
In the case of an infinite uniform Lambertian target of reflectance ρ surf λ and the reflectance of the top of the atmosphere ρ TOA λ at a particular wavelength λ can be written as follows [43]: where ρ atm λ is the atmospheric path reflectance, µ s is the cosine of the view zenith angle (θ s ), µ v is the cosine of the solar zenith angle (θ v ), φ is the relative azimuth angle (the difference between the solar azimuth angle and the satellite azimuth angle), T ↑ λ (µ v ) is the total upward transmission in the direction of the satellite field of view, T ↓ λ (µ s ) is the total downward atmospheric transmission, and S λ is the spherical albedo of the atmosphere (atmospheric backscattering ratio).
On the left-hand side of Equation ( 3), ρ TOA λ can be obtained from satellite observation data.Except for the surface reflectance ρ surf λ , each term on the right-hand side of Equation ( 3) is a function of the aerosol type and the loading optical thickness (τ).Therefore, we can estimate surface reflectance in multiple channels using the three parameters (S λ , ρ atm λ , T ↓ λ (µ s )T ↑ λ (µ v )) and match the AOD value according to the spectral relationship.

Look-Up Table (LUT)
Usually, a set of sun-satellite geometries, atmospheric parameters, and aerosol information is applied to build an LUT, which is then used in the retrieval step to speed up the calculation process.We pre-computed parameters (S λ , ρ atm λ , T ↓ λ (µ s )T ↑ λ (µ v )) using the 6SV radiative transfer model [43,44]  to forward simulation.These parameters were calculated for nine solar zenith angles (from 0 to 80 in increments of 10), nine satellite zenith angles (from 0 to 80 in increments of 10), and 10 relative azimuth angles (from 0 to 180 in increments of 20).In addition, the AOD values were set for 20 aerosol loading, as follows: 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.According to the top of the atmosphere (TOA) reflectance and the atmospheric parameters corresponding to the optical thickness of the continental aerosol and the geometric observation conditions in the LUT, we can calculate the surface reflectance.When comparing the multispectral surface reflectance from the calculation and the surface reflectance from the empirical relationship, and by matching the minimum inversion condition by interpolation, we can retrieve the AOD.

Appropriate Pixel Selection
The DT algorithm is used to retrieve the aerosol properties over cloud-free and dark land surfaces.Therefore, some obstructions such as a clouds, water, snow/ice, and bright pixels need to be identified before AOD retrieval.
We applied a cloud mask method [45] that is based on the spatial variability of reflectance at the TOA to screen out most of the clouds.The visible reflectance test in 0.47 µm can effectively discriminate bright clouds over low reflective surfaces.We can visually interpret the clouds and surfaces [46].Due to basic physical differences between clouds and aerosols, cloud exhibit more complex spatial behavior than aerosols [45].Thus, in the standard deviation of every 3 × 3 pixels for the band 1 (0.46 µm) TOA reflectance (referred to as 3 × 3-STD), clouds show higher values than aerosol.The operational threshold (3 × 3-STD > 0.01, and ρ TOA 0.46 > 0.4) was defined based on the analysis of several cases with different types of clouds and aerosols.Figure 1 shows the true color image (Figure 1a), 3 × 3-STD image in band 1 (Figure 1b) and our cloud mask image (blue color) on 22 September 2016 over the north of China (Figure 1c).corresponding to the optical thickness of the continental aerosol and the geometric observation conditions in the LUT, we can calculate the surface reflectance.When comparing the multispectral surface reflectance from the calculation and the surface reflectance from the empirical relationship, and by matching the minimum inversion condition by interpolation, we can retrieve the AOD.

Appropriate Pixel Selection
The DT algorithm is used to retrieve the aerosol properties over cloud-free and dark land surfaces.Therefore, some obstructions such as a clouds, water, snow/ice, and bright pixels need to be identified before AOD retrieval.
We applied a cloud mask method [45] that is based on the spatial variability of reflectance at the TOA to screen out most of the clouds.The visible reflectance test in 0.47 µm can effectively discriminate bright clouds over low reflective surfaces.We can visually interpret the clouds and surfaces [46].Due to basic physical differences between clouds and aerosols, cloud exhibit more complex spatial behavior than aerosols [45].Thus, in the standard deviation of every 3 × 3 pixels for the band 1 (0.46 µm) TOA reflectance (referred to as 3 × 3-STD), clouds show higher values than aerosol.The operational threshold (3 × 3-STD > 0.01, and ρ 0.46 TOA > 0.4) was defined based on the analysis of several cases with different types of clouds and aerosols.Figure 1 shows the true color image (Figure 1a), 3 × 3-STD image in band 1 (Figure 1b) and our cloud mask image (blue color) on Sep.22, 2016 over the north of China (Figure 1c).In this work, the identification of snow and ice was accomplished by adapting this method using the visible band 2 and near-infrared band 6 [47].The Normalized Difference Snow/ice Index (NDSI) is defined as: NDSI = (ρ 0.51 -ρ 2.30 )/(ρ 0.51 -ρ 2.30 ), (4) where the ρ 0.51 and ρ 2.30 are the AHI TOA reflectances in band 2 and band 6.In this work, the identification of snow and ice was accomplished by adapting this method using the visible band 2 and near-infrared band 6 [47].The Normalized Difference Snow/ice Index (NDSI) is defined as: NDSI = (ρ 0.51 − ρ 2.30 )/(ρ 0.51 − ρ 2.30 ), where the ρ 0.51 and ρ 2.30 are the AHI TOA reflectances in band 2 and band 6.
In addition, the dark pixels are selected based only on their reflectances in near-IR (for MODIS: 2.13 µm) [19,41].To be selected, for AHI, a pixel must fall within the range of ρ 2.30 ≤ 0.25.

Reflectance Relationship Transformation
The signal ρ TOA λ is composed of surface contributions (ρ surf λ ) and atmospheric contributions (ρ atm λ ).We can pre-compute ρ atm λ using the 6SV radiative transfer model.The ground surface reflectances ρ surf λ are difficult to distinguish from the total satellite received signal ρ TOA λ (differentiating the aerosol signal from that of the surface), and the estimation of the surface reflectances is the key factor in aerosol retrieval.The surface reflectance in the visible channel is a function of the surface reflectance in near-infrared channel and is also a function of the scattering angle and the vegetation index [32]: where Θ is the scattering angle, which is defined as follows: where θ s , θ v , and φ are the solar zenith, satellite view zenith, and relative azimuth angles, respectively.In addition, NDVI SWIR is the normalized difference vegetation index in the 1.24 µm channel (MODIS channel 5) and the 2.12 µm channel (channel 7), which are less affected by aerosol.
The function of Equation (3) requires a large number of satellite/ground-based samples to fit.However, as Figure 2 shows, the AERONET sites are located at a few different places, and thus, this distribution is not representative.Correspondingly, the AHI and AERONET data (Table 3) are not sufficient to build the relationship.
The function of Equation ( 3) requires a large number of satellite/ground-based samples to fit.However, as Figure 2 shows, the AERONET sites are located at a few different places, and thus, this distribution is not representative.Correspondingly, the AHI and AERONET data (Table 3) are not sufficient to build the relationship.In this paper, we present a new aerosol retrieval algorithm that applies surface reflectance from MODIS data to retrieve AHI aerosols.Because the differences in their wavelength ranges and spectrum response functions (SRF) (Figure 3) are not great, the surface reflectances can be inherited from MODIS, but they cannot be used directly.In this paper, we present a new aerosol retrieval algorithm that applies surface reflectance from MODIS data to retrieve AHI aerosols.Because the differences in their wavelength ranges and spectrum response functions (SRF) (Figure 3) are not great, the surface reflectances can be inherited from MODIS, but they cannot be used directly.Therefore, we use a transformation method, assuming the reflectances between the MODIS and AHI channels have linear relationship: Therefore, we use a transformation method, assuming the reflectances between the MODIS and AHI channels have linear relationship: where A and B are the coefficients of the function and ρ s_MODIS λ and ρ s_AHI λ are real surface reflectances from MODIS and AHI, respectively, which can be derived by: where Γ(λ) is the SRF (Figure 3); λ 1 and λ 2 are the upper and lower limits of integration, respectively, which is the range of the bands; and, S(λ) is the spectral curve corresponding to central wavelength λ.
In this work, we select 50 kinds of typical features from the spectral data in the ENVI spectrum database (Figure 4) and calculate the reflectance of ground objects in different bands (the blue band, red band, and near-infrared band) of the AHI and MODIS sensors.After differentiation, Equation ( 9) can be rewritten as: Figure 5 shows the linear fitting of reflectance in the three bands, and Table 4 presents the coefficients of A and B, and R 2 .Then, by including reflectance relationship from MODIS [32], we can obtain the AHI surface reflectance relationship: Because the AHI sensor lacks a channel at 1.24 µm, we assume that the value of NDVI SWIR is 0.5 (which is the median of the range in the MODIS algorithm) [32].where A and B are the coefficients of the function and ρ λ s_MODIS and ρ λ s_AHI are real surface reflectances from MODIS and AHI, respectively, which can be derived by: where Γ λ is the SRF (Figure 3); λ 1 and λ 2 are the upper and lower limits of integration, respectively, which is the range of the bands; and, S λ is the spectral curve corresponding to central wavelength λ.
In this work, we select 50 kinds of typical features from the spectral data in the ENVI spectrum database (Figure 4) and calculate the reflectance of ground objects in different bands (the blue band, red band, and near-infrared band) of the AHI and MODIS sensors.After differentiation, Equation ( 9) can be rewritten as: Figure 5 shows the linear fitting of reflectance in the three bands, and Table 4 presents the coefficients of A and B, and R 2 .Then, by including reflectance relationship from MODIS [32], we can obtain the AHI surface reflectance relationship: Because the AHI sensor lacks a channel at 1.24 µm, we assume that the value of NDVISWIR is 0.5 (which is the median of the range in the MODIS algorithm) [32].

Results and Analysis
In this paper, we used the method introduced above to retrieve the AHI AOD product over mainland China.In addition, we extracted AOD_DT, AOD_DB, AOD_V, and AOD_JAXA for comparison and validation.To compare the regional distribution of AOD, AOD_DT, AOD_DB, AOD_NEW and AOD_JAXA are averaged seasonally.In addition, we calculated the difference between AOD_DT, AOD_DB, AOD_NEW, and AOD_JAXA.To quantitatively verify our new algorithm, our AOD_NEW, AOD_JAXA, AOD_DT, AOD_DB and AOD_V were each compared with ground-based observations.In addition, AOD_NEW and AOD_JAXA were also validated at different times from 4 UTC to 9 UTC.Furthermore, we selected data with a one-year span at the XiangHe and Beijing-CAMS AERONET stations for a time series comparison with corresponding satellite data.Finally, AOD_NEW and AOD_JAXA were cross-verified with AOD_DB.

Spatial Distributions of AODs
The spatial and temporal distributions of AOD exhibited great differences.To better illustrate these differences, we calculate the seasonal means and the differences of these data.Figure 6 shows the seasonal mean AOD (computed by the mean of daily data) from the AOD_NEW, AOD_DB and AOD_ DT datasets and their differences.Overall, the spatial distribution of AOD_NEW is more consistent with those of AOD_DB and AOD_ DT.In particular, the difference between AOD_NEW, AOD_DB, and AOD_DT is less than 0.1 in central and eastern China.However, we cannot acquire results over bare land, desert, and other brighter surfaces that do not satisfy the DT condition, such as northwest China.This result illustrates that the new algorithm is more feasible for vegetated cover.As a result of the basic DT algorithm assumption, AOD_NEW is more similar to AOD_DT.For a seasonal comparison, the differences between AOD_NEW and other products are smaller in summer and autumn than in winter and spring.Because of strict cloud mask and pixel selection processing, the AOD_DT results have less effective data under heavy pollution conditions than the AOD_DB product, with similar results for the AOD_NEW product.The deep blue algorithm can perform on days with heavier pollution, accordingly raising the seasonal average and reasonably reflecting the actual atmospheric conditions, which may be why the differences between AOD_NEW and AOD_DB reach 0.4 in parts of the North China Plain (NCP), whereas the differences with AOD_DB are not significant in comparison.Some research [11] has concluded that the applicability of the DB algorithm is better than that of the DT algorithm under pollution conditions.As such, when compared with AOD_DB, our product (with the DT algorithm) has a lower value in all seasons in typical polluted areas in China (such as the Beijing-Tianjin-Hebei area and Chengdu-Chongqing).Figure 7, similar to Figure 6, shows the seasonal mean AOD from the JAXA AHI, MODIS DB, and MODIS DT datasets and their differences.In general, AOD_JAXA exhibits significant underestimation and does not reflect the AOD distribution characteristics.Along the Qinling Mountains-Huaihe River boundary, the value in the north is clearly slightly higher and it is significantly underestimated in the south, where the differences exceed 0.4 in places.Due to the high differences in topography, we can draw the conclusion that surface contributions are the main factor affecting the accuracy of AOD_JAXA.6, shows the seasonal mean AOD from the JAXA AHI, MODIS DB, and MODIS DT datasets and their differences.In general, AOD_JAXA exhibits significant underestimation and does not reflect the AOD distribution characteristics.Along the Qinling Mountains-Huaihe River boundary, the value in the north is clearly slightly higher and it is significantly underestimated in the south, where the differences exceed 0.4 in places.Due to the high differences in topography, we can draw the conclusion that surface contributions are the main factor affecting the accuracy of AOD_JAXA.In Figures 6 and 7, we can see the differences in the seasonal distributions of MODIS/AQUA (DB, DT) AOD, AOD_NEW and AOD_JAXA.In Figure 8, we directly compare AOD_NEW and AOD_JAXA to show the differences more intuitively.The figure shows that the values of AOD_NEW were significantly higher than those of AOD_JAXA, especially in spring and winter in central and eastern China.To a certain extent, in the summer, JAXA products can reflect the distribution of aerosol loading.In Figures 6 and 7, we can see the differences in the seasonal distributions of MODIS/AQUA (DB, DT) AOD, AOD_NEW and AOD_JAXA.In Figure 8, we directly compare AOD_NEW and AOD_JAXA to show the differences more intuitively.The figure shows that the values of AOD_NEW were significantly higher than those of AOD_JAXA, especially in spring and winter in central and eastern China.To a certain extent, in the summer, JAXA products can reflect the distribution of aerosol loading.In Figures 6 and 7, we can see the differences in the seasonal distributions of MODIS/AQUA (DB, DT) AOD, AOD_NEW and AOD_JAXA.In Figure 8, we directly compare AOD_NEW and AOD_JAXA to show the differences more intuitively.The figure shows that the values of AOD_NEW were significantly higher than those of AOD_JAXA, especially in spring and winter in central and eastern China.To a certain extent, in the summer, JAXA products can reflect the distribution of aerosol loading.

The Aerosol Optical Depth Validation Using The AErosol RObotic NETwork (AERONET) Data
In this paper, we select one AERONET site, Beijing at 05:30 (UTC) from July 2016 to June 2017, to calculate the standard deviation of AOD_AERO and AOD_AHI.The correlation between AOD_AERO and AOD_NEW has a relatively high R 2 in Figure 9.The data is dense at low values and it become sparser at a higher optical thickness.As the value becomes larger, the standard deviation becomes larger.In terms of satellite data, this may be due to regional differences.

The Aerosol Optical Depth Validation Using The AErosol RObotic NETwork (AERONET) Data
In this paper, we select one AERONET site, Beijing at 05:30 (UTC) from July 2016 to June 2017, to calculate the standard deviation of AOD_AERO and AOD_AHI.The correlation between AOD_AERO and AOD_NEW has a relatively high R 2 in Figure 9.The data is dense at low values and it become sparser at a higher optical thickness.As the value becomes larger, the standard deviation becomes larger.In terms of satellite data, this may be due to regional differences.

The Overall Comparison
Examining the distribution of AOD in different seasons, we can clearly see that AOD_NEW products, when compared to AOD_JAXA, are closer to the MODIS DT and DB products, which worked successfully for more than a decade.According to above section, we have reason to believe that our results are somewhat more accurate.To better evaluate the accuracy of the product, we need to compare and verify the product the ground-based AOD observations.Because of the temporal resolution (hourly observations for the AHI sensor), the number of samples from MODIS (AQUA) and VIIRS is relatively small.Using the latitude and longitude of the eight AERONET sites, we matched the satellite data with in situ observations from July 2016 to June 2017 in the Chinese region (including Taiwan and Hong Kong).
According to the matching method mentioned in Section 2, the AOD products are matched with ground-based data.Figure 10 shows the comparison and verification results of AOD products, including AOD_NEW, AOD_JAXA, MODIS (DB, DT) AOD, and AOD_V with high quality flag data.The R 2 values were ordered, as follows: AOD_JAXA < AOD_V < AOD_NEW < AOD_DT < AOD_DB.The correlation between AOD_NEW and AOD_AERO has a relatively high R 2 (R 2 = 0.830), and the root-mean-square error (RMSE) is 0.112.However, AOD_NEW is slightly overestimated at low values, which may be due to the aerosol model, which used a single continental aerosol in this paper.The results of AOD_JAXA are the least related to the ground observations and exhibit clear

The Overall Comparison
Examining the distribution of AOD in different seasons, we can clearly see that AOD_NEW products, when compared to AOD_JAXA, are closer to the MODIS DT and DB products, which worked successfully for more than a decade.According to above section, we have reason to believe that our results are somewhat more accurate.To better evaluate the accuracy of the product, we need to compare and verify the product the ground-based AOD observations.Because of the temporal resolution (hourly observations for the AHI sensor), the number of samples from MODIS (AQUA) and VIIRS is relatively small.Using the latitude and longitude of the eight AERONET sites, we matched the satellite data with in situ observations from July 2016 to June 2017 in the Chinese region (including Taiwan and Hong Kong).
According to the matching method mentioned in Section 2, the AOD products are matched with ground-based data.Figure 10 shows the comparison and verification results of AOD products, including AOD_NEW, AOD_JAXA, MODIS (DB, DT) AOD, and AOD_V with high quality flag data.The R 2 values were ordered, as follows: AOD_JAXA < AOD_V < AOD_NEW < AOD_DT < AOD_DB.The correlation between AOD_NEW and AOD_AERO has a relatively high R 2 (R 2 = 0.830), and the root-mean-square error (RMSE) is 0.112.However, AOD_NEW is slightly overestimated at low values, which may be due to the aerosol model, which used a single continental aerosol in this paper.The results of AOD_JAXA are the least related to the ground observations and exhibit clear underestimation relative to AOD_AERO, which is consistent with the comparison of seasonal differences in Figure 8.The comparisons of the AOD_AERO and AOD of the five products show that the AOD from the AHI new algorithm has some reliability.underestimation relative to AOD_AERO, which is consistent with the comparison of seasonal differences in Figure 8.The comparisons of the AOD_AERO and AOD of the five products show that the AOD from the AHI new algorithm has some reliability.

Comparison at Different Times
To verify the stability of AOD_NEW at different times, we extract satellite data and site data from 2:00 to 7:00 (UTC).Figures 11 and 12 shows the comparison scatterplots of AOD_NEW versus AOD_AERO and AOD_JAXA versus AOD_AERO at different times.Overall, the correlation is significantly better in the afternoon than that in the morning for both of the products.AOD_NEW is overestimated, especially between 04:00 to 07:00 (UTC), while AOD_JAXA exhibits significant underestimation.Except for that at 02:00 (UTC), the R 2 of AOD_NEW is higher than that of AOD_JAXA, and we cannot see obvious differences in RMSE.

Comparison at Different Times
To verify the stability of AOD_NEW at different times, we extract satellite data and site data from 2:00 to 7:00 (UTC).Figures 11 and 12 shows the comparison scatterplots of AOD_NEW versus AOD_AERO and AOD_JAXA versus AOD_AERO at different times.Overall, the correlation is significantly better in the afternoon than that in the morning for both of the products.AOD_NEW is overestimated, especially between 04:00 to 07:00 (UTC), while AOD_JAXA exhibits significant underestimation.Except for that at 02:00 (UTC), the R 2 of AOD_NEW is higher than that of AOD_JAXA, and we cannot see obvious differences in RMSE.
underestimation relative to AOD_AERO, which is consistent with the comparison of seasonal differences in Figure 8.The comparisons of the AOD_AERO and AOD of the five products show that the AOD from the AHI new algorithm has some reliability.

Comparison at Different Times
To verify the stability of AOD_NEW at different times, we extract satellite data and site data from 2:00 to 7:00 (UTC).Figures 11 and 12 shows the comparison scatterplots of AOD_NEW versus AOD_AERO and AOD_JAXA versus AOD_AERO at different times.Overall, the correlation is significantly better in the afternoon than that in the morning for both of the products.AOD_NEW is overestimated, especially between 04:00 to 07:00 (UTC), while AOD_JAXA exhibits significant underestimation.Except for that at 02:00 (UTC), the R 2 of AOD_NEW is higher than that of AOD_JAXA, and we cannot see obvious differences in RMSE.

Time Series Comparison of the AErosol RObotic NETwork (AERONET) Data and the Advanced Himawari Imager (AHI) Data
A comparison of different long time series data at certain AERONET stations provides information that can be used to analyze the differences between the AOD products.We select two AERONET sites, Beijing-CAMS (an urban station) and XiangHe (background station), for this comparison.In Figure 13, the three kinds of data (AOD_DB, AOD_AERO, and AOD_NEW) exhibit similar trends, especially at the non-urban XiangHe station.At the Beijing-CAMS station, AOD_NEW also exhibits a uniform trend, but it is somewhat overestimated.As overestimation is observed at the urban site but not at the background station, we can infer that this effect may be caused by surface contributions.Not considering NDVI in the spectral relationship is the main factor.Most of the eight AERONET sites we choose are located in cities, which is why the AOD_NEW in Figure 10 is relatively high when compared to the ground-based sites.The inversion method of the DT algorithm performs better in the countryside than for the brighter surfaces of urban areas, and this result is similar to the current research results for MODIS [48].

Time Series Comparison of the AErosol RObotic NETwork (AERONET) Data and the Advanced Himawari Imager (AHI) Data
A comparison of different long time series data at certain AERONET stations provides information that can be used to analyze the differences between the AOD products.We select two AERONET sites, Beijing-CAMS (an urban station) and XiangHe (background station), for this comparison.In Figure 13, the three kinds of data (AOD_DB, AOD_AERO, and AOD_NEW) exhibit similar trends, especially at the non-urban XiangHe station.At the Beijing-CAMS station, AOD_NEW also exhibits a uniform trend, but it is somewhat overestimated.As overestimation is observed at the urban site but not at the background station, we can infer that this effect may be caused by surface contributions.Not considering NDVI in the spectral relationship is the main factor.Most of the eight AERONET sites we choose are located in cities, which is why the AOD_NEW in Figure 10 is relatively high when compared to the ground-based sites.The inversion method of the DT algorithm performs better in the countryside than for the brighter surfaces of urban areas, and this result is similar to the current research results for MODIS [48].

Time Series Comparison of the AErosol RObotic NETwork (AERONET) Data and the Advanced Himawari Imager (AHI) Data
A comparison of different long time series data at certain AERONET stations provides information that can be used to analyze the differences between the AOD products.We select two AERONET sites, Beijing-CAMS (an urban station) and XiangHe (background station), for this comparison.In Figure 13, the three kinds of data (AOD_DB, AOD_AERO, and AOD_NEW) exhibit similar trends, especially at the non-urban XiangHe station.At the Beijing-CAMS station, AOD_NEW also exhibits a uniform trend, but it is somewhat overestimated.As overestimation is observed at the urban site but not at the background station, we can infer that this effect may be caused by surface contributions.Not considering NDVI in the spectral relationship is the main factor.Most of the eight AERONET sites we choose are located in cities, which is why the AOD_NEW in Figure 10 is relatively high when compared to the ground-based sites.The inversion method of the DT algorithm performs better in the countryside than for the brighter surfaces of urban areas, and this result is similar to the current research results for MODIS [48].

Comparison between Multiple Data Sources
Mature MODIS AOD products, which have been verified around the world and operating robustly for more than 15 years, provide a standard for cross-validating other satellite products.We cross-validated AOD_NEW and AOD_JAXA in three selected regions in China (114.5 14 shows that the correlation between AOD_DB and AOD_NEW is higher than that between AOD_DB and AOD_JAXA.In addition, AOD_NEW is slightly overestimated, especially in low-value regions.AOD_JAXA is relatively low in high-value regions.For some outliers, this result may be caused by the matching between the pixels.

Comparison between Multiple Data Sources
Mature MODIS AOD products, which have been verified around the world and operating robustly for more than 15 years, provide a standard for cross-validating other satellite products.We cross-validated AOD_NEW and AOD_JAXA in three selected regions in China (114.5°E-117.5°E, 38° N-41° N; 112° E-115° E, 31° N-34° N; 113° E-116° E, 23° N-26° N). Figure 14 shows that the correlation between AOD_DB and AOD_NEW is higher than that between AOD_DB and AOD_JAXA.In addition, AOD_NEW is slightly overestimated, especially in low-value regions.AOD_JAXA is relatively low in high-value regions.For some outliers, this result may be caused by the matching between the pixels.

Conclusions
In this work, we introduced a new DT AOD retrieval algorithm for the AHI sensor from July 2016 to June 2017; when compared AOD_NEW with AOD_JAXA, AOD_DT, AOD_DB, and AOD_V over China; and, evaluated the performance of these products with AERONET observations.Most significantly, we used a channel transformation method to revise the surface reflectance relationship from MODIS and adapted this relationship for the AHI sensor.When comparing the spatial distribution of our results with that of MODIS data, we can see that the AHI retrieval with the new algorithm has a better performance than JAXA AOD.Our product has been tested using a one-year time series, and the validation results showed that the product has good correlation and accuracy with AOD_AERO.
Admittedly, our product has some limitations due to simple assumptions in the land cover and the aerosol models.Especially the limitation of the current dark pixel algorithm that is based on the surface reflection relationship can cause systematic error.Although these do not bring very large errors to the existing results, they should be considered as influential factors for aerosol retrieval.They may bring great uncertainty to the results that can be avoided to some extent.
In later research, a dynamic aerosol model can be used to calculate an LUT.In addition, spectral information can be used to determine the best model to improve the AOD retrieval accuracy.Furthermore, a land cover factor (determined by NDVI) needs to be added to the surface reflectance relationship, not being set to a fixed value.As the sensor accumulates more observation data over time, and as more AERONET sites are collected, we can directly obtain the empirical relationship between these two datasets using atmospheric correction.Another solution would be to build a surface reflectance database [10] or a reflectance ratio database [22].
Himawari-8, a new stationary satellite, will be widely used in the observation of aerosol properties in the future.Although little is involved in this study, dynamic information on aerosol loading can be captured due to the high temporal resolution that is offered by this platform.This platform will be highly important in studying aerosol formation and fate, and also have a great role in data assimilation [49,50].In the later work, we will carry out further discussions on high time resolution data based on the results.

Conclusions
In this work, we introduced a new DT AOD retrieval algorithm for the AHI sensor from July 2016 to June 2017; when compared AOD_NEW with AOD_JAXA, AOD_DT, AOD_DB, and AOD_V over China; and, evaluated the performance of these products with AERONET observations.Most significantly, we used a channel transformation method to revise the surface reflectance relationship from MODIS and adapted this relationship for the AHI sensor.When comparing the spatial distribution of our results with that of MODIS data, we can see that the AHI retrieval with the new algorithm has a better performance than JAXA AOD.Our product has been tested using a one-year time series, and the validation results showed that the product has good correlation and accuracy with AOD_AERO.
Admittedly, our product has some limitations due to simple assumptions in the land cover and the aerosol models.Especially the limitation of the current dark pixel algorithm that is based on the surface reflection relationship can cause systematic error.Although these do not bring very large errors to the existing results, they should be considered as influential factors for aerosol retrieval.They may bring great uncertainty to the results that can be avoided to some extent.
In later research, a dynamic aerosol model can be used to calculate an LUT.In addition, spectral information can be used to determine the best model to improve the AOD retrieval accuracy.Furthermore, a land cover factor (determined by NDVI) needs to be added to the surface reflectance relationship, not being set to a fixed value.As the sensor accumulates more observation data over time, and as more AERONET sites are collected, we can directly obtain the empirical relationship between these two datasets using atmospheric correction.Another solution would be to build a surface reflectance database [10] or a reflectance ratio database [22].
Himawari-8, a new stationary satellite, will be widely used in the observation of aerosol properties in the future.Although little is involved in this study, dynamic information on aerosol loading can be captured due to the high temporal resolution that is offered by this platform.This platform will be highly important in studying aerosol formation and fate, and also have a great role in data assimilation [49,50].In the later work, we will carry out further discussions on high time resolution data based on the results.

Figure 2 .
Figure 2. The distribution of AErosol RObotic NETwork (AERONET) stations in China and the number of valid days for the 16 AERONET stations (Table 3) from 7 July 2015 to 1 October 2017.The number of valid days is labeled in different colors.DEM is defined as the Digital Elevation Model.

Figure 2 .
Figure 2. The distribution of AErosol RObotic NETwork (AERONET) stations in China and the number of valid days for the 16 AERONET stations (Table 3) from 7 July 2015 to 1 October 2017.The number of valid days is labeled in different colors.DEM is defined as the Digital Elevation Model.

Figure 4 .
Figure 4. Curve of the spectrum for five kinds of typical features from our 30 selected of features: Brown sandy loam (red line); Lawn grass (blue line); Maple leaves (black line); Aspen leaf (olive line); and, Brown to dark brown loamy sand (orange line).

Figure 4 .
Figure 4. Curve of the spectrum for five kinds of typical features from our 30 selected of features: Brown sandy loam (red line); Lawn grass (blue line); Maple leaves (black line); Aspen leaf (olive line); and, Brown to dark brown loamy sand (orange line).

Figure 4 .
Figure 4. Curve of the spectrum for five kinds of typical features from our 30 selected of features: Brown sandy loam (red line); Lawn grass (blue line); Maple leaves (black line); Aspen leaf (olive line); and, Brown to dark brown loamy sand (orange line).

Figure 5 .
Figure 5. Linear fitting of reflectance in the (a) blue band; (b) red band; and, (c) near-infrared band of AHI and MODIS.

Figure 5 .
Figure 5. Linear fitting of reflectance in the (a) blue band; (b) red band; and, (c) near-infrared band of AHI and MODIS.

19 Figure 6 .
Figure 6.Each column represents a season: from left to right are spring (March, April, and May), summer (June, July, and August), autumn (September, October, and November) and winter (December, January, and February).From top to bottom are the difference between aerosol optical depth (AOD)_NEW and AOD_DB (a-d); the mean value of AOD_DB (e-h); the mean value of AOD_NEW (i-l); the mean value of AOD_DT (m-p); and the difference between AOD_NEW and AOD_DT (q-t).All of the results use full quality data from July 2016 to June 2017.

Figure 6 .
Figure 6.Each column represents a season: from left to right are spring (March, April, and May), summer (June, July, and August), autumn (September, October, and November) and winter (December, January, and February).From top to bottom are the difference between aerosol optical depth (AOD)_NEW and AOD_DB (a-d); the mean value of AOD_DB (e-h); the mean value of AOD_NEW (i-l); the mean value of AOD_DT (m-p); and the difference between AOD_NEW and AOD_DT (q-t).All of the results use full quality data from July 2016 to June 2017.

Figure 7 ,
Figure 7, similar to Figure6, shows the seasonal mean AOD from the JAXA AHI, MODIS DB, and MODIS DT datasets and their differences.In general, AOD_JAXA exhibits significant underestimation and does not reflect the AOD distribution characteristics.Along the Qinling Mountains-Huaihe River boundary, the value in the north is clearly slightly higher and it is significantly underestimated in the south, where the differences exceed 0.4 in places.Due to the high differences in topography, we can draw the conclusion that surface contributions are the main factor affecting the accuracy of AOD_JAXA.

Figure 7 .
Figure 7. Same as Figure 6 but the comparison dataset is AOD_Japan Aerospace Exploration Agency (JAXA).All the results use full quality data from July 2016 to June 2017.

Figure 7 .
Figure 7. Same as Figure 6 but the comparison dataset is AOD_Japan Aerospace Exploration Agency (JAXA).All the results use full quality data from July 2016 to June 2017.

19 Figure 7 .
Figure 7. Same as Figure 6 but the comparison dataset is AOD_Japan Aerospace Exploration Agency (JAXA).All the results use full quality data from July 2016 to June 2017.

Figure 8 .
Figure 8.Each column represents a season, and each row represents (from top to bottom) the following: the difference between AOD_NEW and AOD_JAXA (a-d); the mean value of AOD_NEW (e-h); and the mean value of AOD_JAXA (i-l).All of the results use full quality data from July 2016 to June 2017.

Figure 8 .
Figure 8.Each column represents a season, and each row represents (from top to bottom) the following: the difference between AOD_NEW and AOD_JAXA (a-d); the mean value of AOD_NEW (e-h); and the mean value of AOD_JAXA (i-l).All of the results use full quality data from July 2016 to June 2017.

Figure 9 .
Figure 9. AHI AOD retrieval over land 550 nm as a function of AERONET observations that are collocated in space and time.The value of AOD_AHI: within a 25-km longitude and latitude scope (Beijing: 39.933° N, 116.317°E), more than 25% of the total number pixels are used to calculate the mean.The value of AOD_AERO: using the value of the Beijing site within half an hour before and after 05:30 (UTC) and for values with a corresponding number of ground-based site data points that is no less than 4, and calculating the mean as a transit time site value.The standard deviation in each bin is shown by error bars.

Figure 9 .
Figure 9. AHI AOD retrieval over land 550 nm as a function of AERONET observations that are collocated in space and time.The value of AOD_AHI: within a 25-km longitude and latitude scope (Beijing: 39.933 • N, 116.317 • E), more than 25% of the total number pixels are used to calculate the mean.The value of AOD_AERO: using the value of the Beijing site within half an hour before and after 05:30 (UTC) and for values with a corresponding number of ground-based site data points that is no less than 4, and calculating the mean as a transit time site value.The standard deviation in each bin is shown by error bars.

Table 3 .
The number of valid days from 7 July 2015 to 1 October 2017 for the 32 AERONET stations in China.The eight AERONET stations we used in this work are Beijing, Beijing-CAMS, Chiayi, EPA-NCU, Hong_Kong_PolyU, Lulin, XiangHe and Bamboo.

Table 3 .
The number of valid days from 7 July 2015 to 1 October 2017 for the 32 AERONET stations in China.The eight AERONET stations we used in this work are Beijing, Beijing-CAMS, Chiayi, EPA-NCU, Hong_Kong_PolyU, Lulin, XiangHe and Bamboo.

Table 4 .
Coefficient of slope, intercept and correlation obtained by comparing the MODIS reflectance with the AHI reflectance from the linear fit in Figure5.