Himawari-8-Derived Aerosol Optical Depth Using an Improved Time Series Algorithm Over Eastern China

Abstract

Himawari-8 (H8), as a new generation geostationary meteorological satellite, has great potential for monitoring the spatial–temporal variation of aerosol properties. However, the large amount of spectral data with differing observation geometries require re-formulation of the surface reflectance correction to utilize this new satellite data. This is achieved by using an improved version of the time series (TS) technique proposed by Mei et al., (2012) based on the assumption that the ratio of the surface reflectance in different spectral bands does not change between any two scan times within an hour. In addition, more suitable aerosol models were adopted, based on cluster analysis of local Aerosol Robotic Network (AERONET) data. The improved TS algorithm (ITS) was applied to retrieve the Aerosol Optical Depth (AOD) over eastern China and the results compare favorably with collocated reference AOD data at eleven sun photometer sites (R > 0.8, Root Mean Square Error (RMSE) < 0.2). Comparison with the H8 official AOD product and with MODIS Dark Target (DT)–Deep Blue (DB) combined AOD data shows the good performance of the ITS method for AOD retrieval with different observation angles.

1. Introduction

Aerosols are important for climate and climate change, due to their reflection and absorption of solar radiation and their effects on cloud properties [1]. They are also important for human health, especially for people with respiratory problems and lung diseases [2,3]. Furthermore, they play a role in atmospheric chemistry, while at the same time their variability is very large, which severely affects cloud microphysics [4,5]. Hence, it is important to obtain accurate information on the occurrence of aerosols and their properties. One way to obtain this information on large spatial scales is the use of satellites, using sensors with an initiative-designed task for aerosol retrieval. Many sensors designed for land and sea surface temperatures are also suitable for aerosol retrieval, such as the European Space Agency’s (ESA) along-track scanning radiometers (ATSR-2 and AATSR) [6,7,8].

Sensors suitable for aerosol retrievals are Moderate Resolution Imaging Spectroradiometer (MODIS) [9], Visible infrared Imaging Radiometer (VIIRS) [10], Polarization and Directionality of the Earth’s Reflectances (POLDER) [11] and the Chinese directional polarimetric camera (DPC) [12], with different characteristics for aerosol retrieval, such as wavebands at wavelengths ranging from the ultra-violet (UV) to the thermal infrared (TIR; MODIS, VIIRS), two or more viewing angles (ATSR, POLDER), and polarization (POLDER, DPC). All of these sensors are attached to sun-synchronous, polar-orbiting satellites, which implies that the revisit time for each point on Earth is one day or more (depending on swath width). Therefore, at best one or two observations per day are available [13,14,15]. In this study, we focus on the use of geostationary satellites, and in particular the advanced Himawari imager (AHI) onboard the Himawari-8 (H8) satellite. Geostationary satellites observe only a certain part of the earth, however at a temporal frequency of several observations each hour during the daytime. This provides information on the diurnal variation of aerosol properties, which is important for air quality studies and health effects, as well as to study aerosol transport. In the current study, we only consider aerosol retrieval over land.

Only a few algorithms have reached a state of maturity, rendering them suitable for the operational production of aerosol data and the creation of climate data records [16,17]. The most crucial issue for the retrieval of aerosol properties using passive sensors is the effective separation of the contributions from the underlying surface and atmospheric components to the reflectance measured at the top of the atmosphere (TOA), especially because the combination is nonlinear and under-constrained [18,19].

The wide variety of land surface types results in variability of the surface appearance from dark forest to bright desert, snow, or ice covered surfaces, with an associated variation of surface reflectance [20]. Due to the huge spatial–temporal variability in aerosol concentration and the complex mixing of aerosol types, aerosol retrieval is difficult [4,21,22]. Therefore, to retrieve aerosol information from satellite observations, effective separation of the surface reflectance and atmospheric path radiance contributions to the reflectance measured at the TOA is necessary. Furthermore, satellite observations do not provide enough information for a detailed characterization of the aerosol properties, and therefore assumptions of the aerosol type are required.

To effectively correct for the surface contributions to the TOA reflectance, several methods have been proposed using different assumptions. For sensors with a single view, the accuracy of the aerosol signal separated from the TOA signal is largely based on an accurate determination of the surface reflection [19,23]. For MODIS, the dark target (DT) method has been developed [24], which uses the real-time relationship between the reflectance in the visible and shortwave infrared bands over dark dense vegetation. However, the DT method applies only to a limited range of surface albedos [9], while for more reflecting surfaces other methods are needed. This idea was also later used in the retrieval of sensor Advanced Baseline Imager (ABI) (on American stationary satellite Geostationary Operational Environmental Satellites - R Series (GOES-R)) [25]. The deep blue (DB) method [26] was developed for use over brighter surfaces, which utilizes wavelengths in the UV spectral region where the surface appears dark. These methods have also been applied to other sensors, such as VIIRS, Sea-Viewing Wide Field-of-View Sensor (SeaWiFS), and Advanced Very High Resolution Radiometer (AVHRR) [27]. DB and DT data products are, in principle, complementary, covering different types of surfaces, while MODIS Collection 6 and 6.1 data product include a combined DT and DB AOD product [28,29,30]. Other methods use physical models to represent the angular variation of the surface reflectance, as described by the bidirectional reflectance distribution function (BRDF). A widely used BRDF data set is available from MODIS observations using the Ross–Li model based on images measured over 16 consecutive days [31]. For geostationary satellites, it is relatively easier to obtain enough observations. Govaerts et al. [32,33,34] used the Rahman–Pinty–Verstraete (RPV) model to simultaneously retrieve the average AOD and surface BRDF reflectance of Meteosat Second Generation–Spinning Enhanced Visible and Infra-red Imager (MSG–SEVIRI) within a day. Because of the high temporal resolution of geostationary satellites, their data can facilitate surface and aerosol research.

In the current study, we focus on the use of the H8–AHI geostationary satellite. For aerosol retrieval, Ge et al. [35] used the relationship between the reflection in the visible and mid-infrared bands of H8–AHI, following the MODIS and ABI algorithms [24,25]. However, due to the long observation period of a geostationary satellite, this method has some shortcomings in establishing a strong relationship between the surface reflectance in the visible spectral band and the TOA reflectance in the infrared band. She et al. [36] used the MODIS BRDF dataset [31] to jointly retrieve the hourly AOD and the surface reflectance, but this method also suffers from some errors due to fixing the shape of the surface BRDF. In the H8–AHI official aerosol retrieval algorithm (further referred to as the “official algorithm”), a predefined minimum surface method was used from a historical dataset with hourly equivalent BRDF information [37] (similar to Geostationary Ocean Color Imager (GOCI) [38,39]). This method may introduce some uncertainties by synthesizing the optimal albedo from historical data, resulting in obvious deviations in some areas compared with AOD measurements from AERONET sites [40,41].

In addition to surface assumptions, there are still some difficulties in deriving AOD from atmospheric signals. The aerosol microphysical parameters determine the scattering and absorption of light by particles, so the use of incorrect aerosol parameters directly leads to incorrect AOD retrievals. In general, aerosol properties can be obtained from the analysis of long-term measurements [21,32]. However, for simplicity, many algorithms use existing aerosol models and do not take into account the spatiotemporal variations of actual aerosol types [35,42].

Mei et al. [43] presented a time series (TS) method for MSG–SEVIRI data, which is based on the assumption that the surface albedo in one pixel remains unchanged during three consecutive scans. With predefined fixed aerosol types, this method obtained the AOD by utilizing a numerical minimization routine. The algorithm makes full use of the characteristics of high temporal resolution and does not need to obtain accurate surface reflectance; therefore, it can be applied to various surface types. In addition, the algorithm can also obtain the optimal aerosol type from predefined aerosol types by limiting the spatial variation of aerosol types.

In theory, the TS method is especially suitable for geostationary satellites. However, the original TS algorithm is not effective on H8–AHI due to method accuracy and inappropriate aerosol parameters [44,45,46]. Therefore, this study attempts to improve the TS algorithm in the following five aspects: (1) re-derive the basic transfer function and consider gas correction and Rayleigh scattering coupling; (2) modify the cost function to better consider the ratio of surface changes; (3) extend the observation period of the TS method to retrieve more retrievals; (4) set up aerosol models through cluster analysis to get more accurate results; (5) use optimal estimates for fast iteration with an efficient constraint. The results obtained using the improved time series (ITS) algorithm have a cross-comparison with ground-based data, the MODIS AOD product, and the H8–AHI official AOD dataset. A detailed description of the retrieval methodology is provided in Section 2, where the data are also introduced. The comparison between retrieval results is discussed in Section 3. Conclusions are presented in Section 4.

2. Retrieval Strategy and Data

2.1. Advanced Himawari Imager

The advanced Himawari imager (AHI) is a 16-channel, multispectral imager onboard the geostationary satellite Himawari-8 (H8), developed by the Japan Meteorological Agency (JMA). H8–AHI perform full-disk observations over East Asia and the western Pacific at 10 min intervals. The first 6 bands (0.47 μm–2.2 μm) from the visible to the mid-infrared can be used to capture the radiance of atmospheric aerosols and the earth’s surface at a maximum spatial resolution of 1 kilometer in East Asia at 10 min intervals. To facilitate processing and utilization, JAXA released real-time L1B spectral data with 5 km resolution in NC format [15,47].

In the current study on the development of the ITS algorithm, we use the 5 km resolution L1B spectral dataset available from the JAXA website (http://www.eorc.jaxa.jp/ptree/index.html). Real-time L2 cloud products with the same spatial–temporal resolution are used to select cloud-free pixels [48,49]. L1B TOA albedo values at 470, 510, 640, 870, and 2250 nm measured at 10 min intervals from 09:00 to 16:00 China Standard Time (CST, i.e., UTC + 8:00) over the study area are used as input. The study period is from August 2018 to May 2019.

Complementary to satellite observations, local AOD measurements are made using ground-based sun photometers. The data from sun photometer networks, such as the global Aerosol Robotic Network (AERONET) [50,51] or the Chinese Sun–Sky Radiometer Observation Network (SONET) [52], are more accurate than satellite products, and therefore are used as a reference for the validation of satellite-retrieved AOD [35,36]. In the current study, ground-based AOD measurements from eleven AERONET and SONET sun photometer sites (see

Table 1. Locations of the eleven sun photometer sites used in the current study, including the surface type and period for which data were available.

Sites	Longitude (°)	Latitude (°)	Altitude (m)	Period	Surface
Beijing-CAMS	116.38	39.98	92	2018.8–2019.05	Urban
Beijing	116.38	40.01	59	2018.8–2019.05	Urban
Xianghe	116.96	39.75	36	2018.8–2019.05	Rural
Xuzhou	117.14	34.22	60	2018.8–2019.05	Suburb
Taihu	120.22	31.42	20	2018.8–2018.10	Urban
Yanqihu	116.67	40.41	100	2018.8–2019.05	Rural
Jiaozuo	113.25	35.19	113	2018.8–2019.05	Urban
Songshan	113.10	34.54	475	2018.8–2019.05	Woodland
Hefei	117.16	31.91	36	2018.8–2019.05	Suburb
Nanjing	118.96	32.12	52	2018.8–2019.05	Suburbs
Shanghai	121.48	31.28	85	2018.8–2019.05	Urban

for details) are used for validation. The four AERONET sites are Beijing-CAMS, Xianghe, Xuzhou, and Taihu. The seven SONET sites are Beijing, Yanqihu, Jiaozuo, Songshan, Hefei, Nanjing, and Shanghai. AERONET provides aerosol data at three quality levels: level 1.0 (unscreened), level 1.5 (cloud-screened and quality controlled), and level 2.0 (quality-assured) (http://aeronet.gsfc.nasa.gov). The quality of SONET AOD data corresponds to that of AERONET Level 1.5 [45,46]. Because L2.0 AERONET data for 2018 were not yet available at the time of the current study, Level 1.5 AERONET V3 AOD data were used together with corresponding SONET data.

For comparison, the H8–AHI official version 2.1 level 2.0 5km AOD product [53] and the MODIS C6.1 10 km AOD data (DT and DB combined dataset, DT–DB) from Terra and Aqua were selected. The MODIS 10 km product is widely used and can be obtained from the LAADS website (http://ladsweb.nascom.nasa.gov/data). MODIS is in a sun-synchronous polar orbit, with equator crossing times at about 10:30 (Terra satellite) and 13:30 (Aqua satellite) China Standard Time (CST).

2.2. Basic Method

The ITS algorithm needs to construct the cost function based on the surface reflectance. The ITS method uses the same two-stream approximate radiative transfer model as the original TS method (classic formula derived by Liou et al.) [54]. This model decomposes the pure atmospheric scattering processes in spherical coordinates into upward and downward flux streams ($F^{\uparrow },F^{\downarrow }$; $\tau$ represent the atmospheric optical thickness at any position):

$$\begin{array}{l}\frac{d{F}^{\uparrow }(\tau )}{d\tau }={\gamma }_1{F}^{\uparrow }(\tau )-{\gamma }_2{F}^{\downarrow }(\tau )-{\gamma }_3\varpi F_0\exp (-\frac{\tau }{{\mu }_0})\\ \frac{d{F}^{\downarrow }(\tau )}{d\tau }={\gamma }_2{F}^{\uparrow }(\tau )-{\gamma }_1{F}^{\downarrow }(\tau )+(1-{\gamma }_3)\varpi F_0\exp (-\frac{\tau }{{\mu }_0})\end{array}$$

(1)

In Equation (1), $F_0$ and ${\mu }_0$ represent the solar flux and cosine of its direction angle, while the coefficients ${\gamma }_1,{\gamma }_2,{\gamma }_3$ are only related to the single scattering albedo (ssa, $\varpi$) and the asymmetry factor (asy, $g$) of the atmosphere layer. To integrate these formulas, the boundary conditions need to be clear. Then, considering the surface albedo and the albedo at TOA, $A,A^{\prime }$, the boundary conditions for upward and downward fluxes at the top and bottom of the atmosphere are reset as Equation (2).

$$\frac{F\uparrow \left(\tau =0\right)}{{\mu }_0{F}_0+F\downarrow \left(\tau =0\right)}=A^{\prime }$$

(2a)

$$F\downarrow \left(\tau =0\right)=0$$

(2b)

$$\frac{F\uparrow \left(\tau ={\tau }_0\right)}{F\downarrow \left(\tau ={\tau }_0\right)}=A$$

(2c)

In Equation (2), $\tau$ is the atmospheric optical depth, where $F$ (arrow up and down with $\tau$ = 0) indicates the fluxes at TOA, and likewise $F$ (arrow up and down with $\tau$ = ${\tau }_0$) indicates the fluxes at the surface. Based on known boundary conditions and atmospheric parameters, the surface albedo $A$ can be derived from TOA reflection $A^{\prime }$. After mathematical deformation and simplification, the final formula can be formulated as Equation (3):

$$A=\frac{\left(a-c{A}^{\prime }\right)e^{k\tau }+\left(b+c\Gamma A^{\prime }\right)\Gamma e^{-k\tau }+\left({\Gamma }^2-1\right)G^{+}{e}^{-\frac{\tau }{{\mu }_0}}}{\left(a-c{A}^{\prime }\right)\Gamma e^{k\tau }+\left(b+c\Gamma A^{\prime }\right)e^{-k\tau }+\left({\Gamma }^2-1\right)\left(G^{-}+c\right)e^{-\frac{\tau }{{\mu }_0}}}$$

(3)

In Equation (3), since no new variables were introduced, $a,b,c,\Gamma$ are derived from the parameters in Equation (1), so they are only related to atmospheric parameters and incident solar radiation. They are functions of the incoming solar radiation flux ($F_0$), solar zenith angle (${\mu }_0$), single scattering albedo (ssa,$\tilde{\omega }$), and asymmetry factor (asy, $g$) [55,56]. All the parameters are used by Mei et al., (2012), which are the deformations of the two-stream parameters (a detailed explanation of the symbols and the equation derivation are presented in the Appendix A).

The parameters of solar flux are easily obtained, so this formula indicates that for known atmospheric parameters (aerosol and gas), the surface albedo can be calculated from the TOA reflectance measured by H8–AHI.

For a clear sky, we assume that the atmospheric optical depth ($\tau$) consists of two parts (Equation (4a)): the molecular Rayleigh scattering (${\tau }_R$) and the aerosol optical depth (${\tau }_A$). An approximate expression for the Rayleigh scattering, Equation (4b), is sufficiently accurate for most applications in remote sensing [1]. The aerosol optical depth (AOD) varies with wavelength, as described by the Ångström turbidity Formula (4c) [1]. Over a small area and within a short period time, the aerosol properties (represented by the wavelength exponent $\alpha$) can be considered to be invariant, while the concentration of the aerosol particles (given by the Ångström turbidity coefficient $\beta$) may change.

$$\tau ={\tau }_A+{\tau }_R$$

(4a)

$${\tau }_R^{\lambda }(\infty )=0.00879{\lambda }^{-4.09}$$

(4b)

$${\tau }_A=\beta {\lambda }^{-\alpha }$$

(4c)

Due to the mixing of all atmospheric particles, the optical properties of the atmospheric layer use Equations (5a) and (5b) to describe the atmospheric coupling single scattering albedo (${\tilde{\omega }}^{\prime }$) and coupling asymmetric factor ($g^{\prime }$) by optical depth weighting [57,58].

$${\tilde{\omega }}^{\prime }=({\tilde{\omega }}_A·{\tau }_A+{\tau }_R)/({\tau }_A+{\tau }_R)$$

(5a)

$$g^{\prime }=(g_A·{\tau }_A·{\tilde{\omega }}_A)/({\tau }_A·{\tilde{\omega }}_A+{\tau }_R)$$

(5b)

Additionally, gas absorption should be considered. Similar to the procedure adopted by MODIS, we use the National Centers for Environmental Prediction (NCEP) data to correct for the effect of gas absorption (H₂O, CO₂, and O₃) on the measured TOA reflectance before retrieval is attempted [59]. The NCEP data include the 1° × 1° global meteorological analysis data. For each pixel in the L1B data, the corrected reflectance $A^{\prime }$ is given by Equation (6), which has a scale factor $T_{\lambda }^{gas}$ based on wavelength, air mass factor, and weighting coefficients (Appendix B and Appendix C; benefiting from the similar band between AHI and MODIS, here we directly use MODIS parameter values).

$${A^{\prime }}_{\lambda }=T_{\lambda }^{\mathrm{gas}}{A}_{\lambda }^{L1B}$$

(6)

In the case of strongly absorbing aerosols and an optically thin layer, a negative surface albedo may result from the numerical correction due to the high forward scatter peaks of atmospheric particles (an inherent defect of the two-stream approximation). Therefore, we also adjust this case by using a $\delta$ function that removes the amount of scattering of the forward peak while ensuring the consistency of parameters $\tau$, $\varpi$, and $g$ [54].

2.3. Aerosol Types

Aerosol properties are important and inappropriate aerosol type can result in inaccurate retrieval results. In the literature, different models have been proposed. For example, Govaerts et al., (2010) suggested a model consisting of six aerosol types (three for spherical and three for non-spherical particles) for use in MSG–SEVIRI AOD retrieval [32]. This model was also used by Mei et al., (2012). Derived from AERONET global data, the aerosol model used in the MODIS DT algorithm is another good model [28]. However, considering the scope of the study area, the optical properties of the global model do not agree well with local areas (Figure 1a,b). The scatterplots of the ssa versus asy at an AERONET wavelength of 675 nm, together with the values for Govaerts and MODIS aerosol models, show that these models are not completely representative for the actual aerosol properties over eastern China, especially for the absorbing spherical particles (ssa is on the edge of the cloud of data points or beyond).

Therefore, it is more suitable to construct aerosol types from long-term local measurements and accumulated meteorological data. However, a suitable aerosol model needs to consider the microphysical properties of the particles and local climatology, which is very complicated. For the two-stream approximation method in this study, the optical absorption parameter (ssa) and size parameter (asy) derived from the aerosol microphysical properties are sufficiently accurate to define different aerosol types [22,43]. Therefore, a simple k-means clustering method can be used for aerosol classification based on ssa and asy for AERONET observations in different bands. The theoretical understanding of this statistical method is incomplete; however, it can avoid explicit assumptions regarding aerosol physics in other aerosol models.

To this end, AERONET V3 level 2.0 aerosol data (best quality) at 440, 500, 675, and 870 nm spectral bands from more than 10 local AERONET sites over eastern China were collected since 2010. For ssa and asy in different bands, the k-means clustering method only needs to specify the number of types to divide the different types. Considering the average difference of each type and other models, 5 types are finally determined. For comparison, the results at 675 nm are presented in Figure 1c, with each of the five clusters coded in different colors. The values for other H8–AHI wavelengths were obtained by interpolation using quadratic fitting, which are shown in Figure 1d–f and

Table 2. Aerosol types resulting from the cluster analysis described in the text: single scattering albedo values (ssa, $\tilde{\omega }$) and asymmetry parameters (asy, $g$) at four different wavelengths.

Model $\tilde{\mathit{\omega }}/\mathit{g}$	470 nm	510 nm	640 nm	870 nm
1	0.941/0.743	0.946/0.736	0.963/0.711	0.962/0.696
2	0.839/0.697	0.83/0.688	0.814/0.664	0.785/0.659
3	0.944/0.70	0.946/0.689	0.953/0.653	0.947/0.632
4	0.89/0.704	0.891/0.696	0.895/0.672	0.88/0.66
5	0.895/0.673	0.897/0.66	0.904/0.618	0.889/0.60

.

In addition, AOD conversion between bands for different types can use the Ångström exponent (AE). AE is provided directly in the AERONET dataset and the highest frequency AE value in the histogram of each cluster is determined as the final AE. Using this AE value, the AOD at different wavelengths can be obtained using Equation (4c). Then, Equation (3) becomes a function that only contains AOD as unknown.

2.4. The Core Strategy

In the original TS algorithm presented by Mei et al. [1], the surface reflectance was assumed to be invariant between three consecutive observations within 30 min. There are three main disadvantages: (1) If the surface reflectance values in different bands have multiple differences, then the square sum will further expand this difference, resulting in a band with high reflectance becoming dominant in the cost function. For example, the 640 nm band of the dark pixel in the original DT algorithm is twice as large as the 470 nm band [24]. (2) The assumption of invariant surfaces leads to an uncertain time period for maintaining this invariance. In fact, although the surface reflectance changes slowly with time, it still cannot be considered constant due to the complex BRDF [33,34]. (3) Due to the cloud mask, the requirement for effective observations in three consecutive times also greatly reduces the retrieved pixels. Therefore, the cost function needs to be reconstructed.

For multi-viewing-angle instruments, the observations at different angles allow for simultaneous retrieval of surface reflectance and AOD [60]. The ratio of the directional reflectance values for different viewing directions, here referred to as the K-ratio, is approximately constant at different wavelengths [1,19,61]. This principle can also be applied to geostationary satellites based on the reflectance variation caused by different angles of the sunlight over a period of time. An effective cost function (Equation (7)) is established, from which the most suitable result should satisfy the minimum sum of differences between the K-ratios in different bands.

The cost function of ITS can be evaluated using the following equation:

$${\epsilon }_{final}=\underset{m=5}{\min }\left\{{\displaystyle \sum _{i=1}^n{\displaystyle \sum _t{\left(A_{{\lambda }_{t1},i}^m/A_{{\lambda }_{t2},i}^m-A_{{\lambda }_{t1},i+1}^m/A_{{\lambda }_{t2},i+1}^m\right)}^2}}\right\}$$

(7)

where $m$ is the aerosol type (the total number of aerosol types is 5) defined for observations at time $t1,t2$ (which represents any scans within one hour) and different wavebands ($i$ represents the spectral band, $n$ is the total number of bands equal to 4). The different K-ratios can be calculated for different AOD combinations to determine the minimum value for the cost function $\epsilon$. In this cost function, the subtractions of different K-ratios take on an order of wavelength. In fact, the pair can be arbitrarily chosen here and the result is not much different (because the values of K-ratios change slowly within an hour). By changing aerosol type, five values $\epsilon$ can be obtained, with the smallest one corresponding to the most probable AOD and aerosol type in the retrieval area.

In addition, by using K-ratios invariance, the observation intervals can also be expanded. The H8–AHI official document (constructed surface reflectance data set for each hour during daytime using the second minimum method) and She et al. also recommend one hour [36]. Considering the spatial–temporal variation of aerosol properties, in this paper, if there are two effective scans in an hour, Equation (6) can be used. Therefore, in the absence of sufficient observations due to the real weather, we can also get enough data. Then, an optimal constraint estimation is used to search for the best result within the bounds of real aerosol and surface conditions (assuming that the surface reflectance of the 640 nm visible band is less than the TOA reflectance of the 2250 nm middle infrared band).

The AERONET wavelengths used to define the aerosol types in this study (Figure 1 and Table 1) are close to the H8–AHI visible and NIR bands, so the optical properties ssa and asy of H8–AHI spectra (470, 510, 640, and 860 nm) can be obtained by quadratic fitting from AERONET data. For a single pixel, the selection of aerosol types may have some uncertainties due to the model accuracy, the aerosol type error, and other factors. An important constraint for the final aerosol type is used: the aerosol properties do not change during one hour and over a 1° × 1° spatial grid [22,43]. With this scheme, the aerosol type for pixels within an area can be updated by counting the aerosol type that appears most.

2.5. Execution Steps

Figure 2 shows the flow chart of the ITS aerosol retrieval algorithm. For each series of cloud-free pixels from more than two H8–AHI scans within an hour, the preliminary result consists of an aerosol type, AODs, and surface albedos (for each scan), which are retrieved in steps 1 and 2. For the pixel series in each 1° × 1° grid, the most frequent aerosol type for about 400 pixels would be set as the final aerosol type for that grid (step 3). That aerosol type is used for the final processing of the 1° × 1° image (step 4) and provides the final result (step 5). This procedure imposes a spatial constraint on the aerosol types, which to some extent ensures the accuracy of the final result.

2.6. Study Area

The study area is located between 105° to 125°E and 25° to 45°N, which covers Eastern China, the Bohai Sea, and part of the Huanghai Sea (Figure 3). Many megacities are located in this area, meaning it is strongly influenced by anthropogenic atmospheric pollution [2,5]. Because of the differences between land and ocean algorithms, for this study, we only consider the retrieval of aerosols over land and validate the results by comparison with data from eleven ground-based AOD reference sites in Eastern China. The locations of these AERONET and SONET sites are marked in Figure 3 in yellow and red, respectively.

2.7. Evaluation Metrics

Spatial–temporal matching criteria are the most important criteria for the comparison of different datasets. The MODIS dataset was resampled to the same resolution as H8–AHI (5 km) by using the nearest neighbor interpolation method. However, the AODs from the MODIS product and the ITS-derived AODs are retrieved at 550 nm, the H8–AHI official product is retrieved at 519 nm, and both are different from the wavelength of the ground-based reference data. Therefore, all AOD data were converted to 0.55 µm using Equation (3c). For the convenience of comparison, the overpass times of MODIS Terra and Aqua were defined to 10:30 CST and 13:30 CST (slightly different from the actual local time).

The AOD retrieved for one pixel using the ITS method is only related to the reflectance of that same pixel in two consecutive moments, which is independent of other surrounding pixels. In order to verify the collocations from different datasets accurately, a spatial–temporal matching method used the average values of ground-based measurements during a certain time period (10 min, ± 5 min around the satellite overpass time) and the average of the satellite-retrieved AOD over an area of 5 × 5 pixels around the sun photometer site (for stability, more than 10 valid values in 25 pixels were included the calculation) was used [10,43].

For the evaluation of the H8–AHI ITS algorithm, scatterplots were made of the retrieved AOD versus reference AOD data. Four parameters were used as metrics to evaluate the performance of the ITS method: (i) the correlation coefficient (R); (ii) the percentage of ITS retrievals falling within the expected error (EE) ($\pm (0.05+0.15\tau )$), which applies to MODIS data over land and has been adopted by others [9,10]; (iii) mean difference ($MD={\displaystyle \sum ({\tau }_{retrieved}-{\tau }_{groundbased})/n}$ describes the mean bias between two datasets; (iv) root mean square error ($RMSE=\sqrt{{\displaystyle \sum {({\tau }_{retrieved}-{\tau }_{groundbased})}^2/n}}$) refers to the satellite product’s absolute deviation from ground-based AOD).

3. Results and Analysis

3.1. Validation against AERONET Measurements

Scatterplots of satellite-retrieved AOD products versus collocated ground-based AOD data are presented in Figure 4. For comparison, the H8–AHI L2 official AOD product is shown in the left column, the MODIS product is shown in the middle column, and the ITS method is shown in the right column. Overall, the three datasets compare favorably with the ground-based reference measurements. However, at most sites, the data retrieved using the ITS algorithm have better correlation and lower RMSE than the official H8–AHI product, and at all sites more collocations fall within the EE than the official H8–AHI product. The results from the ITS algorithm are overall not as good as those from MODIS, but at some sites (Nanjing, Taihu, and Shanghai) the ITS statistics are better than those from MODIS.

The official algorithm has taken into account the hourly BRDF characteristics, but the result is the worst in the three datasets with small correlation coefficients (often <0.8) and large RMSEs (>0.32). Especially at Yanqihu, Beijing, Beijing-CAMS, and Xianghe sites, the official product at low AOD (AERONET AOD<0.15) is much too high. In these cases, the results from the ITS method are much better than those from the official product, with a smaller number of overestimations, a larger R (>0.8), and a smaller RMSE (<0.20). Because both algorithms use the same cloud screening product, cloud contamination cannot be the reason for the higher AOD of the official product. Levy et al., (2010) concluded that surface assumptions tend to dominate AOD retrieval when there is low aerosol loading (AOD<0.15) [9]. Therefore, we conclude that the surface albedo assumption with the constant K-ratio used in the ITS algorithm provides a better surface correction than the minimum reflectance method used in the official product. However, as mentioned in Table 1, the underlying surface of Beijing and Beijing-CAMS belongs to the urban surface, whereas the Yanqihu and Xianghe sites are located in rural areas.

For Beijing and Beijing-CAMS, the ITS method partially overestimates the AOD (which is not encountered in other sites with an urban underlying surface, such as Jiaozuo and Shanghai), indicating that the surface correction in the ITS method may need further consideration. At the Xuzhou site, the distributions of three datasets are similar and the correlation coefficient of H8–AHI official product is larger than that for the ITS and MODIS. However, the collocations from ITS have the smallest RMSE and more collocations fall in the EE lines. Except for the two sites in Beijing, the ITS-derived data from other sites have a similar distribution (the upper or lower deviations of the EE envelope) from the MODIS product. When the underlying surface is a different type, there is no significant difference in the results, which indicates that the ITS algorithm can successfully retrieve AOD irrespective of surface type and the aerosol model is representative of the actual situation. H8–AHI provides much more data than MODIS due to the high temporal frequency, which fully reflects the advantages of the ITS algorithm applied to geostationary satellites. At the Taihu site, due to the weather, the number of successful retrievals from the ITS algorithm is less than that from the official product because two validated cloud-free pixels are needed in one hour. This results in the rejection of more pixels than the H8–AHI official method, which uses a single cloud-free observation.

3.2. Time Series Analysis

The frequent observations of geostationary satellites allow for monitoring of the diurnal evolution of the AOD. Figure 5 shows the AOD time series over a period of several days for selected sites, and for one site (Taihu) a full-time series of about 3 months. To this end, collocations between three satellite datasets and ground-based observations were selected from Figure 4, for which large numbers of AOD retrievals with different aerosol loading rates are available. The data in Figure 5 show that for different aerosol loadings, the ITS-derived AOD generally traces the ground-based observations well. Absolute differences are less than 0.2 and vary throughout the day. This indicates that the algorithm can capture small aerosol changes that vary with solar irradiation and the different underlying surfaces.

At the Beijing-CAMS (panel a) and Xianghe (panel b) sites, the performances of the H8–AHI official product and ITS method are different from those shown in Figure 4. The official product has a large overestimation under low aerosol loading, which is also reflected in Figure 5. Although the reference data show an almost constant AOD, those from the official product increase substantially with time and remain high. This indicates that the surface hypothesis and aerosol type used in the official algorithm do not always provide the correct result. In contrast, the ITS-derived AODs at the Beijing site are often closer to the reference, although high AOD is also observed in that data. At the Xianghe site, the ITS-derived AOD traces the reference data quite well and deviations are small, regardless of the level of aerosol loading. This illustrates the applicability of the ITS algorithm for multiple surface types. However, the performance in Beijing indicates that the ITS algorithm still needs to be further processed under some complex surfaces.

The Jiaozuo site (line c) is located in another urban area. Compared with Beijing, the AODs derived by the ITS agree well with the ground-based observations, and the daily variation is very smooth. This indicates that the ITS algorithm has a good performance, which it also has over some urban surfaces.

For the Taihu site, a full time series of about 3 months (4 August to 10 October 2018) is presented. The performance of the H8–AHI official algorithm is good at the Taihu site, with more successful retrievals than the ITS method. The ITS-derived AOD and the official product both trace the reference data well, often with a small underestimation by the ITS and a larger overestimation by the official product.

3.3. Spatial–Temporal Distributions of ITS-Derived Data

Maps of the AOD over eastern China, derived using the ITS method, are shown in Figure 6 for each hour from 09:00 to 16:00 LST on October 04, 2018. The maps show the spatial variation over the study area and the evolution during the day, with a strong increase shown over almost the whole study area, and in particular over the North China Plain. This fully demonstrates the advantages of geostationary satellites for aerosol research.

3.4. Comparison with Official H8–AHI and MODIS AOD Products

This section focuses on the spatial distribution of the ITS-retrieved AOD and its performance compared with that from the H8–AHI official product and the MODIS products. AOD maps for each of these three products, retrieved from observations on October 5, 2018, during the time of the Terra and Aqua overpasses (about 10:30 and 13.30 CST, respectively), are presented in Figure 7. Seasonal mean results including September and October 2018 are compared in Figure 8. The spatial distributions of the three datasets in Figure 7 show similar patterns but the values are very different. At 10:30 CST, the coverage of the official product (panel b) is very sparse. In contrast, in the afternoon, the H8–AHI official product (panel e) AOD is much larger than that from the ITS or MODIS. ITS-derived AOD (c and f) is much closer to the MODIS AOD (a and d), although somewhat lower. The aggregated data of the three dataset the MODIS Aqua overpass time (about 13:30 CST) can further exhibit the difference. To calculate the means, only the pixels for which the number of retrievals is larger than 2 (at least two cloud-free days) were used to calculate the monthly means. The ITS and MODIS AOD mean maps show a similar spatial pattern, while the AOD of the official product is overall much higher. This observation is consistent with the scatterplots in Figure 4. Comparison of the ITS-derived AOD maps with those from MODIS shows a qualitative similar spatial distribution, however the ITS-derived AOD is somewhat lower and patchier. This may be due to the internal accuracy difference of the two-stream model in the ITS algorithm.

For further comparison of the ITS-derived AOD and the MODIS product, the pixel-by-pixel differences between the ITS and Terra (about 10:30 CST) products are plotted in Figure 9a, while the difference plot for ITS Aqua (about 13:30 CST) is presented in Figure 9b. Figure 9a shows that in the northern part of the study area, the ITS-derived AOD is generally higher than the MODIS AOD, whereas in the southern part of the area the ITS-derived AOD is smaller. Some areas are slightly overvalued (red), but undervalued areas (blue) are more common. The coverage of positive and negative biases does not have a clear range boundary or trend between morning and afternoon. The scatterplots in Figure 9c shows that the bulk of the data is centered around the identity line, indicating good quantitative agreement between ITS and MODIS Terra for AOD up to about 0.5. However, large deviations occur, especially overestimation by a factor of 2–3. The comparison with MODIS Aqua is much lower; especially for low AOD, the ITS-derived results have an obvious underestimation and there seems to be little relation between the ITS and MODIS Aqua AOD during the afternoon overpass. The difference in the performance of the ITS algorithm between the morning and afternoon indicates that further improvements need to be made.

3.5. Aerosol Type Analysis

The ITS method can simultaneously retrieve AOD and aerosol-type products (from five predefined types). Most aerosol retrievals focus on AOD accuracy [10,62], but aerosol-type products are also a research topic [60,63,64]. In this section, the preliminary results from the ITS method are discussed.

With the two-stream approximate method in the ITS algorithm, different aerosol types can be retrieved for adjacent pixels with similar aerosol loading (accuracy limitations). This is physically not acceptable, and therefore for each 1° × 1° grid cell the aerosol type, which occurs in the largest number of pixels, is selected as the most probable type for that grid cell. It appears that most of the time and in most areas, moderately absorption (type 3; Figure 1 and Table 1) occurs most frequently. This type also occurs most frequently according to AERONET statistics. This may be due to the relatively small study area where anthropogenic aerosols dominate.

Figure 10 shows the time series of three aerosol parameters (AOD, ssa, and asy at 670 nm) derived from AERONET measurements and retrieved by the ITS algorithm for single pixels at the Beijing site. The quality of the ITS dataset at this site is not as good as other sites (Figure 4), which is representative of studying the relationship between aerosol types and AOD. Overall, the ITS-derived AODs have a qualitatively similar variation to the AERONET data. The most corresponding points, where ssa is equal to 0.953 and asy is equal to 0.653, belong to the third type. The largest differences between the ITS and AERONET AOD occur when the ssa and asy are also different from the reference data, such as on January 9 and February 28. The AOD variation is similar to that of the reference AERONET data, however the AOD values diverge when the ssa deviation is small and asy deviation is large, such as the date before January 17 and the period between February and March. When ssa and asy are in good agreement with the actual aerosol properties from ground-based measurements, this also applies to ITS-derived AOD (after March).

As expected, the aerosol type directly affects the accuracy of the AOD retrievals. As shown in Figure 10 (top panel), the ITS-derived AODs do not compare well with the reference data until February, whereas after February the comparison is good and the ITS AOD traces the AERONET AOD. The overestimation in the winter may be due to the bright surface at that time of the year, together with low AOD, as observed by AERONET [9]. During episodes with higher AOD, the ITS data compare favorably with the AERONET data.

4. Discussion and Conclusion

H8–AHI, which has an unprecedented high temporal resolution, has the capability to capture the dynamic variations of aerosols, and thus provide an important contribution to atmospheric environmental monitoring [65,66]. The use of H8–AHI data for aerosol retrieval requires an effective separation of the effects of the surface and the atmosphere, with the reflectance measured by the sensor at the top of the atmosphere. To this end, the use of a fast and accurate algorithm derived by Mei et al., (2012) was proposed and further improvements were made for use with H8–AHI data. Compared with the official algorithm and the MODIS DT algorithm, this algorithm can make full use of the characteristics of geostationary satellites. The generic nature of this algorithm renders it suitable for real-time operational use. Simultaneously retrieved aerosol types give the algorithm wider application prospects. Validation with independent reference data from AERONET and comparison with MODIS AOD show the good performance at the ITS-retrieved AOD. The ITS algorithm performance is better than that of the current official H8–AHI aerosol product (smaller RMSE and more points within the EE at all sites) and close to the MODIS product. The aerosol models are used in the algorithms, which were developed for this area using the clustering method for AERONET data in the study area. The evaluation of the retrieved aerosol types shows a good agreement between the retrieval aerosol type and the real aerosol type.

However, improvements to this new algorithm are still needed. Comparison with AERONET measurements at the Beijing-CAMS and Beijing sites shows an obvious overestimation at low aerosol loading (AOD<0.2), although the ITS results are better than the official product. Additionally, the ITS performs better at other urban sites than at the Beijing sites. This may be explained from the residual contribution of the surface albedo and the differences in retrieval algorithms between ground-based instruments and satellite sensors (which may be insensitive to very thin aerosol layers [54]). The overestimation of the AOD for low aerosol loading at the Beijing sites has also been reported in other studies [35,36] and requires further analysis. In addition, because this algorithm uses only two effective scans in an hour, it lacks a strict surface constraint over one or more days [67].

The use of accurate aerosol types is crucial. The aerosol type retrieved by the ITS algorithm has a good correlation with the real situation, and the influence of the absorption as expressed by the ssa and the asy on AOD was found (Figure 10). However, this algorithm only considers these two parameters, so it should be improved in the future to cover a variety of atmospheric properties for a specific area, such as eastern China [62]. On the other hand, the degrees of freedom available from satellite observations are insufficient to provide detailed information on the aerosol properties, which may lead to a simple formulation of aerosol models, providing a priori estimates, which are refined in the retrieval processing using the observations at different wavelengths [33,34].

The uncertainty factors mentioned above will be addressed in a future study to improve the performance of the ITS algorithm for atmospheric monitoring using the H8–AHI satellite.