SAHARA: A Simplified AtmospHeric Correction AlgoRithm for Chinese gAofen Data: 1. Aerosol Algorithm

Abstract

The recently launched Chinese GaoFen-4 (GF4) satellite provides valuable information to obtain geophysical parameters describing conditions in the atmosphere and at the Earth’s surface. The surface reflectance is an important parameter for the estimation of other remote sensing parameters linked to the eco-environment, atmosphere environment and energy balance. One of the key issues to achieve atmospheric corrected surface reflectance is to precisely retrieve the aerosol optical properties, especially Aerosol Optical Depth (AOD). The retrieval of AOD and corresponding atmospheric correction procedure normally use the full radiative transfer calculation or Look-Up-Table (LUT) methods, which is very time-consuming. In this paper, a Simplified AtmospHeric correction AlgoRithm for gAofen data (SAHARA) is presented for the retrieval of AOD and corresponding atmospheric correction procedure. This paper is the first part of the algorithm, which describes the aerosol retrieval algorithm. In order to achieve high-accuracy analytical form for both LUT and surface parameterization, the MODIS Dark-Target (DT) aerosol types and Deep Blue (DB) similar surface parameterization have been proposed for GF4 data. Limited Gaofen observations (i.e., all that were available) have been tested and validated. The retrieval results agree quite well with MODIS Collection 6.0 aerosol product, with a correlation coefficient of R² = 0.72. The comparison between GF4 derived AOD and Aerosol Robotic Network (AERONET) observations has a correlation coefficient of R² = 0.86. The algorithm, after comprehensive validation, can be used as an operational running algorithm for creating aerosol product from the Chinese GF4 satellite.

1. Introduction

The retrieval of aerosol optical properties is a difficult task due to the coupling of surface and aerosol information observed by the satellite Top of Atmosphere (TOA) reflectance. Currently, the determination of aerosol types, including size distribution, refractive index, and aerosol profile, complicates the problem, especially for the aerosol types with strong absorption properties [1,2]. The cloud twilight zone effect [3] is also a challenging task in order to obtain accurate aerosol properties [4]. An aerosol retrieval algorithm is instrument-dependent, especially for the estimation of surface properties. Currently, the majority of aerosol retrieval algorithms like Moderate Resolution Imaging Spectroradiometer (MODIS) Dark-Target (DT), MODIS Deep Blue (DB) are implemented utilizing the Look-Up-Table (LUT) approach [5,6]. The LUT method speeds up the retrieval procedure compared to full radiative transfer calculations. However, accurate analytical/empirical equations provide a unique ability to undertake fast aerosol retrievals, which are especially useful for high-spatial resolution aerosol retrievals [7,8].

Several attempts have been proposed in past decades for fast aerosol retrievals with analytical equations. The MultiAngle Implementation of Atmospheric Correction (MAIAC) algorithm is a semi-analytical algorithm using Green function by utilizing the analytical calculation of the surface bidirectional reflectance factor (BRF), which is proven to be an effective way for the operational aerosol algorithm [9]. Katsev et al. [10] and Seidel et al. [11] used analytical solutions of radiative transfer theory to speed up the AOD retrievals, and the retrievals using MEdium Resolution Imaging Spectrometer (MERIS) instrument show very promising results. Another way is to utilize the Kaufman et al. [12] equation which was firstly proposed by Chandrasekhar [13] and tries to parameterize different terms in the radiative transfer equation, including atmosphere path reflectance, transmittance and spherical albedo. The newly development of Simplified high resolution MODIS Aerosol Retrieval Algorithm (SARA) [7] utilized the parameterization with the aerosol properties derived from a local urban Aerosol Robotic Network (AERONET) station and the surface reflectance from MOD09GA level-2 daily surface reflectance product. The algorithm has been validated and can be further used for the high spatial resolution AOD retrieval with enough local a-priori information.

Besides the parameterization of the radiative transfer process, another essential issue is how to describe different aerosol physical/chemical properties, which is directly linked to the aerosol type in aerosol remote sensing. The accuracy of aerosol type selection has significant influence in atmospheric correction and surface reflectance retrieval. Previous studies show that different aerosol model assumptions may cause approximately 10% differences in the estimated spectral surface reflectance [1,2]. A popular method is to derive regional-representative aerosol information using local in-situ measurements like AERONET. From the global scale, the main aerosol types include biomass burning aerosol over Amazon and South Africa, urban/industry aerosol types over North Africa, Europe and Eastern Asia, dust aerosol types over Africa and maritime aerosol over ocean regions [14,15,16]. The investigation shows that single scattering albedo is a key parameter to describe the positive/negative radiative forcing of atmospheric aerosol while asymmetry factor and AOD shows the magnitude of these effects [16,17]. Thus, we can try to parameterize the aerosol micro-physical properties using those key parameters, and further establish functions to describe the properties of particle size distribution, shape and components [17,18,19]. The selection of optimal aerosol properties which provide good accuracy and high-speed, is another important issue for fast and accurate aerosol retrieval [20]. The MODIS DT algorithm provides the seasonal/geo-location dependent aerosol type information [5]. One of the advantages is that this treatment tries to link all aerosol physical/chemical information to the main retrieval parameter, AOD, and avoid introducing new un-known parameters. Radiative calculations (atmosphere path reflectance, transmittance and spherical albedo) for three pre-defined MODIS aerosol types can be further parameterized for the fast and accurate aerosol retrieval.

Most satellite aerosol products have a spatial resolution of 10 km, which was initially designed for climate applications [5]. However, as the applications like small-scale aerosol events observations using aerosol products have been highlighted in recent years, providing aerosol products with higher spatial resolution becomes a very important issue [21]. The standard aerosol product like MODIS DT provides an operational aerosol product with a spatial resolution of 3 km besides the standard 10 km product in MODIS Collection 6.0 (C6) [21,22]. The MAIAC algorithm provides AOD product at a spatial resolution of 1 km [9]. Many other studies also focus on obtaining even higher spatial resolution AOD data product (with spatial resolution higher than 1 km). Zhang et al. [23] presented an aerosol retrieval algorithm for Geostationary Ocean Color Imager (GOCI) instrument with spatial resolution of 500 meter based on the assumption that change of AOD is the main factor for change of satellite-observed TOA compared to surface reflectance within a short period. Li et al. [24] used a synergetic method by combing Chinese Huanjing-1 (HJ-1) CCD and MODIS data to retrieve AOD with spatial resolution of 100 meter. Zhang et al. [25] utilized the pre-defined relationship between surface reflectance at 0.47 µm and TOA Reflectance of 1.6 µm channels of HJ-1 data to obtain high-resolution AOD (100 meter). Vermote et al. [26] proposed a “ratio map” method to retrieve AOD with spatial resolution of 30 meter for Landsat 8 Operational Land Imager (OLI) instrument. The “ratio map” is a 5.5 km auxiliary database for the surface reflectance ratio between channel 1/2 and channel 4 of the OLI instrument calculated by MODIS surface reflectance product. Other attempts have been made to achieve high spatial resolution AOD from Landsat data [27,28,29] beside the classical method as mentioned in [26]. Some other recent progress of European aerosol research focusing on AATSR dual-view aerosol retrieval can be found in the Climate Change Initiative—Aerosol project [30]. The inversion from AERONET for aerosol and surface reflectance inversion were reported in the AERONET Inversion Products document and can also be found in [31].

Chinese high-resolution earth observing system (CHEOS) has planned a series of high spatial resolution satellite observations named GaoFen (GF) series. GaoFen-4 (GF4), launched on 29 December 2015, which will be used in this paper, is the first geostationary satellite in the GF series. There have been limited attempts at the retrieval of aerosol properties using GF series. Wang et al. [32] applied the modified DB algorithm to GF1 Wide Field of View (WFV) data to retrieve high spatial resolution AOD over Eastern China. Bao et al. [33] proposed an AOD retrieval algorithm for GF1 data based on the modified dynamic relationship of the red/blue surface reflectance. To our best knowledge, there has been no publication related to the retrieval of aerosol properties using GF4 data till now.

In this paper, a Simplified AtmospHeric correction AlgoRithm for Chinese gAofen data (SAHARA) is presented. The main idea is to propose parameterizations for the radiative transfer equation [12] based on the MODIS DT aerosol types [5]. A similar idea to DB has been implemented for the surface parameterization. A statistical relationship between GF4 visible channels (0.47 µm and 0.67 µm) and Rayleigh-corrected TOA reflectance at 0.87 µm for different surface types over China has been obtained using MODIS surface reflectance data product. Then an iteration procedure using Levenberg–Marquardt algorithm [34] is used to retrieve the optimal AOD.

The paper has been organized in the following structure. The details of spectral, spatial and temporal characteristics of GF4 instrument have been listed in Section 2. The surface reflectance and aerosol type parameterization, aerosol retrieval algorithm are included in Section 3. The retrieval results and corresponding analysis/validation are in Section 4.

2. GaoFen-4 (GF4) Data

CHEOS intends to construct an independent and advanced earth observation system providing satellite observations with high spatial, temporal and high spectral resolution [35]. GF4 has six channels between visible and mid-infrared band listed as

Table 1. Instrument characteristics of GaoFen-4 (GF4) camera (provided by China’s center for resources satellite data and application) and Moderate Resolution Imaging Spectroradiometer (MODIS) instrument.

GF4 Band Number	GF4 Band-Width (µm)	Resolution (m)	MODIS Band Number	MODIS Band-Width (µm)	Resolution (m)
1	0.45–0.9	50 m
2	0.45–0.52	50 m	4	0.459–0.479	500 m
3	0.52–0.60	50 m	3	0.545–0.565	500 m
4	0.63–0.69	50 m	1	0.620–0.670	250 m
5	0.76–0.90	50 m	2	0.841–0.876	250 m
6	3.50–4.10	400 m

. It has a spatial resolution of 50 meter for visible channels and 400 meter for infrared channel, with a swath width of 400 km. GF4 satellite is centered at 105.6°N and can observe China and the surrounding areas by point control. The high temporal/spatial resolution enables GF4 to be widely used in resources, environment, agriculture, disaster prevention and reduction. Table 1 shows the instrument characteristics of GF4. GF4 instrument has similar bands as MODIS observations. The difference of band width between GF4 and MODIS can be found. However, the analysis shows very good agreement between TOA reflectance between GF4 and MODIS for selected wavelengths for SAHARA. The correlations coefficients and slopes for selected wavelengths used for SAHARA are above 0.9, indicating the possibility to use MODIS product to derive prior-knowledge for SAHARA.

3. Method

In this section, key steps including surface reflectance and aerosol type parameterization are discussed in detail. The MODIS DT cloud screening algorithm without cirrus cloud test (1.38 µm test) has been used with the current retrieval together and a visual-check according to the RGB images. The surface parameterization is a MODIS DB similar approach [6] while the aerosol types are parameterized by adapting the MODIS DT algorithm [5].

3.1. Surface Parameterization

Surface parameterization is one of the key issues for aerosol retrieval. In the SAHARA algorithm, we follow a similar idea to the MODIS DB aerosol retrieval algorithm [6]. Statistical relationships of surface reflectance between different wavelengths were obtained. To achieve accurate estimation of the underlying surface reflectance, we performed a detailed statistical analysis for each surface type according to the MODIS Land cover product (MCD12C1) with spatial resolution of approximately 5 km [36], in which land surface is classified into 16 classes. The land cover map for the region of (70°–140°E, 15°–55°N) is depicted in Figure 1. The details of surface reflectance estimation for each surface type are described below.

Previous attempts show the possibilities of using MODIS data product for surface parameterization and atmospheric parameters retrieval for instruments like Landsat (similar to GF4) [26,37]. One year (2008) cloud-free MODIS C6 TOA reflectance and corresponding MODIS surface reflectance product (MOD09A1) were collected over the study region as shown in Figure 1. Previous studies show that certain relationships between the surface reflectance at visible channels (0.67 and 0.47 µm) and reflectance at near-infrared/ infrared channels (e.g., 0.86, 1.6, 2.1 and 3.75 µm) can be found for certain surface types depending on instruments [5,6,38,39,40]. Figure 2 shows an example of the relationships between the spectral surface reflectance at visible and near-infrared channels for different surface types over China. According to Figure 2, a linear relationship between the spectral surface reflectance at the visible and near-infrared channels for most surface types can be found, which agrees with previous studies listed above. The surface properties are affected by vegetation amount or greenness [41] illustrated by NDVI’ (defined as Equation (1)). The relationships between visible and near-infrared channels is also partly determined by NDVI’.

NDVI’ = (RCR_0.86 − RCR_0.67)/(RCR_0.86 + RCR_0.67)

(1)

where RCR is Rayleigh Corrected TOA Reflectance for given wavelengths.

In order to avoid potential aerosol contamination in the analysis, match-ups with AODs at 0.55 µm were constrained to values smaller than 0.1 (obtained from MODIS C6 aerosol product) and the scattering angle between 110° and 160° were used in the statistical analysis [38]. The contribution of Rayleigh scattering for channel 0.86 µm was removed from the TOA reflectance following DB algorithm.

Figure 3 shows the linear relationships between collocated MODIS surface reflectance in visible channels and corresponding TOA reflectance data (after removing the Rayleigh scattering contribution) in 0.86 µm channel over China during autumn 2008. Figure 3 confirms that this relationship is a function of NDVI’, and the slopes of surface reflectance from 0.86 to 0.67 μm are more sensitive to the NDVI’ values than the 0.86 to 0.47 μm slopes (clearly decrease of slopes with increase of NDVI for slope_0.67/0.86 compared to slope_0.47/0.86). Besides NDVI’, Levy et al. [41] and Wu et al. [42] found that the relationship for visible and near-infrared also depends on scattering angle. In the statistical analysis, the surface reflectances at the visible channels were determined by Rayleigh-corrected TOA reflectances at 0.86 µm (RCR_0.86), NDVI’ and scattering angle (Ψ) for each surface type.

Based on the analysis above, we established an approach to estimate the surface reflectance at visible channels (0.47 μm and 0.67 μm) determined by the Rayleigh-corrected TOA reflectances at the near-infrared channel (0.86 μm), NDVI’ and scattering angle values using the following formulas:

R_0.67 = slope_0.67/0.86 × RCR_0.86 + yint_0.67/0.86

(2)

slope_0.67/0.86 = a₀ × NDVI’ + a₁×Ψ + a₂

(3)

yint_0.67/0.86 = b₀ × Ψ + b₁

(4)

R_0.47 = slope_0.47/0.86 × RCR_0.86 + yint_0.47/0.86

(5)

slope_0.47/0.86 = c₀ × NDVI’ + c₁ × Ψ + c₂

(6)

yint_0.47/0.86 = d₀ × Ψ + d₁

(7)

where the R_0.67 and R_0.47 are the estimated reflectance, slope_0.67/0.86 and slope_0.47/0.86 are the slope of the regressed relationships of 0.67 vs. 0.86 and 0.47 vs. 0.86, yint_0.67/0.86 and yint_0.47/0.86 are the offset of the regression relationships. a₀, a₁, a₂, b₀, b₁, c₀, c₁, c₂, d₀, and d₁ are the regression coefficients determined by least squares fitting to the real surface reflectance.

For each surface type defined in Figure 1, the regression coefficients were derived seasonally for three NDVI’ groups similar to the MODIS DB algorithm: 0.10 < NDVI’ ≤ 0.20; 0.20 < NDVI’ ≤ 0.50; 0.50 < NDVI’, respectively. The coefficients for each land type and NDVI’ group for four seasons are displayed in

Table 2. Surface Reflectance Coefficients for 0.47 µm for March-April-May (MAM).

Land Cover	0.10 < NDVI’ < 0.20	0.20 ≤ NDVI’ < 0.50	0.5 ≤ NDVI’
1	0.456, −2.450, −0.321, 1.730, −0.168	−0.053, −0.743, −0.337, 1.329, 0.020	−0.186, −0.358, −0.747, 1.912, 0.089
2	−0.951, −1.494, −1.250, 3.953, 0.414	0.170, −1.195, 0.563, 0.085, −0.130	0.354, −0.414, 1.063, −1.610, −0.189
3
4
5	−0.099, −1.886, 0.061, 1.358, −0.011	0.026, −0.432, 0.090, 0.181, −0.010	0.060, −0.499, 0.192, 0.065, −0.031
6		0.273, −1.445, 0.423, −0.499, −0.094
7	−0.006, −0.890, −0.072, 0.562, 0.021	0.583, −0.762, 0.684, −1.429, −0.202
8	−0.596, −0.573, −0.793, 3.101, 0.201	0.100, −0.377, 0.281, −0.254, −0.045	−0.010, −0.363, 0.022, 0.319, 0.003
9
10	0.029, −1.089, 0.059, 0.446, −0.013	0.028, −0.575, 0.091, 0.315, −0.017	0.012, −0.271, 0.052, 0.198, −0.007
11
12	−0.115, −1.219, −0.062, 1.208, 0.010	0.006, −0.500, 0.066, 0.308, −0.005	0.005, −0.433, 0.083, 0.236, −0.004
13	−0.046, −2.090, 0.001, 1.282, −0.013	0.113, −1.019, 0.152, 0.429, −0.050
14	−0.128, −1.445, 0.216, 1.163, −0.045	−0.039, −0.592, −0.057, 0.636, 0.012
15
16	−0.015, −2.396, −0.037, 0.806, 0.014	0.080, −1.502, −0.397, 1.542, −0.015

,

Table 3. Surface Reflectance Coefficients for 0.47 µm for June-July-August (JJA).

Land Cover	0.10 < NDVI’ < 0.20	0.20 ≤ NDVI’ < 0.50	0.5 ≤ NDVI’
1			−0.115, −0.546, −0.211, 1.023, 0.047
2		0.008, −1.519, 0.187, 0.679, −0.002	0.168, −0.678, 0.314, −0.091, −0.076
3
4			0.167, −0.554, 0.522, −0.733, −0.075
5		0.094, −1.448, 0.268, 0.638, −0.067	−0.028, −0.571, −0.024, 0.605, 0.010
6
7	0.201, −1.350, 0.370, −0.433, −0.063
8		0.011, −1.650, 0.202, 1.010, −0.052	−0.124, −0.789, −0.079, 0.963, 0.034
9
10	−0.071, −1.848, −0.123, 1.098, 0.023	−0.039, −0.554, −0.054, 0.651, 0.013	0.021, −0.403, 0.058, 0.251, −0.009
11
12	−0.260, −2.270, −0.232, 1.886, 0.096	0.226, −1.249, 0.444, −0.062, −0.109	−0.189, −0.652, −0.296, 1.338, 0.073
13		−0.104, −1.362, −0.200, 1.471, 0.037
14		−0.563, −1.524, −0.674, 2.946, 0.205	−0.618, −0.703, −0.789, 2.491, 0.262
15
16	0.069, −1.567, 0.133, 0.270, −0.021	−0.009, −0.466, −0.005, 0.403, 0.015

,

Table 4. Surface Reflectance Coefficients for 0.47 µm for September-October-November (SON).

Land Cover	0.10 < NDVI’ < 0.20	0.20 ≤ NDVI’ < 0.50	0.5 ≤ NDVI’
1	0.362, −2.306, 1.662, −1.471, −0.244	0.009, −0.890, 0.106, 0.424, −0.004	0.001, −0.344, 0.032, 0.287, 0.000
2		−0.095, −1.425, −0.025, 1.187, 0.025	−0.004, −0.365, 0.009, 0.343, 0.002
3
4		−0.012, −0.439, −0.009, 0.455, 0.005	−0.227, −0.466, −0.573, 1.775, 0.100
5	0.010, −1.577, 0.291, 0.816, −0.066	−0.012, −0.446, 0.007, 0.406, 0.006	−0.011, −0.418, 0.001, 0.402, 0.006
6		−0.041, −0.501, −0.151, 0.795, 0.021
7	0.028, −0.469, 0.051, 0.221, 0.005	−0.061, −0.488, −0.124, 0.679, 0.041
8	−0.134, −2.022, 0.041, 1.201, 0.057	0.020, −1.193, 0.178, 0.609, −0.028	−0.016, −0.404, −0.004, 0.394, 0.008
9
10	0.010, −0.934, −0.017, 0.505, 0.002	−0.008, −0.595, 0.005, 0.470, 0.005	−0.003, −0.368, 0.029, 0.287, 0.003
11	−0.031, −2.031, 0.562, 0.527, −0.030	0.018, −0.993, 0.031, 0.565, −0.011
12	−0.048, −1.480, 0.044, 0.847, 0.003	−0.021, −0.501, −0.013, 0.473, 0.011	−0.029, −0.503, −0.026, 0.529, 0.012
13	−0.046, −1.490, 0.044, 0.891, 0.011	−0.269, −0.580, −0.692, 2.104, 0.122
14	−0.384, −1.422, 0.044, 1.144, 0.153	−0.035, −0.550, −0.031, 0.591, 0.011	−0.010, −0.422, 0.004, 0.380, 0.006
15
16	0.023, −1.344, 0.056, 0.413, 0.003	−0.034, −0.364, −0.046, 0.444, 0.033

,

Table 5. Surface Reflectance Coefficients for 0.47 µm for December-January-February (DJF).

Land Cover	0.10 < NDVI’ < 0.20	0.20 ≤ NDVI’ < 0.50	0.5 ≤ NDVI’
1	−0.125, −1.788, −0.217, 1.542, 0.045	0.010, −0.959, 0.145, 0.448, −0.010	0.010, −0.365, 0.049, 0.258, −0.004
2		0.068, −0.847, 0.225, 0.237, −0.040	0.008, −0.334, 0.033, 0.265, −0.004
3
4		0.047, −0.550, 0.254, 0.010, −0.027
5	−0.003, −2.118, 0.407, 0.434, −0.026	−0.004, −0.760, 0.077, 0.563, −0.008	−0.009, −0.386, 0.007, 0.349, 0.006
6		−0.015, −0.557, 0.013, 0.438, 0.009
7	0.059, −0.874, 0.114, 0.200, −0.015	−0.040, −0.610, −0.029, 0.427, 0.044
8	−0.104, −1.499, 0.081, 1.068, 0.018	−0.003, −0.619, 0.027, 0.453, 0.000	0.004, −0.382, 0.033, 0.294, −0.001
9
10	0.004, −0.836, 0.090, 0.420, −0.007	−0.014, −0.560, 0.014, 0.471, 0.006	0.018, −0.317, 0.087, 0.132, −0.007
11	−0.061, −2.128, −0.019, 1.184, 0.010	−0.074, −0.696, −0.097, 0.768, 0.036
12	−0.108, −2.426, 0.059, 1.090, 0.025	0.011, −0.609, 0.094, 0.311, −0.006	−0.009, −0.365, 0.015, 0.353, 0.002
13	−0.071, −2.326, −0.018, 1.125, 0.023	−0.037, −0.941, 0.013, 0.694, 0.015	−0.012, −0.432, 0.001, 0.427, 0.005
14		−0.018, −0.572, −0.005, 0.508, 0.007	−0.007, −0.425, 0.003, 0.406, 0.002
15
16		−0.079, −0.365, −0.134, 0.603, 0.063

,

Table 6. Surface Reflectance Coefficients for 0.67 µm for MAM.

Land Cover	0.10 < NDVI’ < 0.20	0.20 ≤ NDVI’ < 0.50	0.5 ≤ NDVI’
1	0.0036, −1.657, −0.116, 1.353, 0.010	0.027, −1.551, −0.033 1.124, 0.003	0.063, −0.542, 0.196, 0.049, −0.0088
2	−0.159, −3.625, −0.121, 1.765, −0.079	−0.251, −0.69, −0.443, 1.617, 0.129	−0.019, −0.661, −0.127, 0.902, 0.028
3
4		−0.281, −1.230, −0.302, 1.791, 0.113
5	−0.081, −1.512, −0.196, 1.496, 0.044	−0.161, −1.071, −0.336, 1.679, 0.082	−0.066, −0.631, −0.205, 1.033, 0.049
6	0.1435, −4.187, 0.358, 0.624, −0.055
7	0.0289, −1.909, −0.046, 1.294, −0.005	−0.081, −0.788, −0.236, 1.550, 0.035
8	−0.0364, −1.963, −0.151,1.500, 0.027	−0.257, −1.287, −0.467, 2.030, 0.123	−0.123, −0.544, −0.351, 1.215, 0.080
9
10	0.0217, −1.2854, 0.069, 0.934, −0.008	−0.0002, −1.075, −0.0536, 1.136, 0.005	−0.043, −0.708, −0.217, 1.111, 0.038
11
12	−0.0228, −2.366, −0.083, 1.399, 0.020	−0.091, −1.244, −0.234, 1.529, 0.049	−0.080, −0.662, −0.218, 1.151, 0.046
13	−0.124, −0.240, −0.329, 1.571, 0.067	−0.0525, −1.233, −0.216, 1.456, 0.034
14	0.008, −2.698, −0.0189, 1.466, −0.013	−0.172, −0.976, −0.392, 1.793, 0.082	−0.0115, −0.60, −0.119, 0.831, 0.023
15	−0.0645, −1.530, −0.191, 1.563, 0.037	0.051, −1.400, −0.0214, 1.058, −0.006	0.032, −0.520, 0.065, 0.435, −0.002
16	0.020, −1.657, −0.116, 1.353, 0.010	−0.001, −1.349, −0.084, 1.234, 0.0136	0.059, −0.846, −0.041, 0.841, −0.009

,

Table 7. Surface Reflectance Coefficients for 0.67 µm for JJA.

Land Cover	0.10 < NDVI’ < 0.20	0.20 ≤ NDVI’ < 0.50	0.5 ≤ NDVI’
1		−0.027, −0.378, −0.371, 1.141, 0.037	0.032, −0.553, 0.006, 0.442, 0.007
2		0.061, −1.330, 0.024, 1.023, −0.024	−0.063, −0.545, −0.196, 0.910, 0.050
3
4			−0.242, −0.459, −0.425, 1.405, 0.121
5		0.091, −1.237, 0.061, 0.956, −0.049	−0.051, −0.409, −0.130, 0.644, 0.041
6			0.144, −0.737, 0.254, 0.055, −0.049
7	−0.077, −1.976, −0.189, 1.725, 0.033	−0.029, −1.351, −0.143, 1.362,0.026
8		−0.041, −1.451, −0.151, 1.480,0.013	−0.120, −0.404, −0.268, 0.954, 0.071
9
10	0.0494, −1.524, 0.017, 1.076, −0.013	0.011, −1.344, −0.060, 1.161, 0.008	−0.033, −0.636, −0.131, 0.870, 0.034
11		0.155, −1.337, 0.464, 0.114, −0.075	−0.112, −0.600, −0.067, 0.668, 0.064
12	−0.359, −2.430, −0.441, 2.403, 0.145	0.135, −0.974, 0.230, 0.257, −0.044	−0.013, −0.647, −0.080, 0.734, 0.025
13		0.016, −1.118, −0.070, 1.015, 0.010	−0.056, −0.687, −0.151, 0.952, 0.043
14		0.045, −1.479, 0.009, 1.018, −0.008	−0.096, −0.490, −0.207, 0.865, 0.064
15
16	0.027, −1.585, −0.032, 1.197, −0.002	0.005, −1.434, −0.083, 1.254, 0.011	0.102, −0.859, 0.086, 0.586, −0.028

,

Table 8. Surface Reflectance Coefficients for 0.67 µm for SON.

Land Cover	0.10 < NDVI’ < 0.20	0.20 ≤ NDVI’ < 0.50	0.5 ≤ NDVI’
1	0.156, −2.234, 0.369, 0.429, −0.078	0.003, −0.860, −0.115, 0.963, 0.014	0.009, −0.433, −0.182, 0.794, 0.018
2		−0.0622, −1.245, −0.149, 1.378, 0.021	−0.032, −0.419, −0.128, 0.636, 0.040
3
4		−0.098, −0.616, −0.366, 1.450, 0.056	−0.039, −0.654, −0.248, 1.105, 0.038
5	0.005, −1.455, 0.147, 1.089, −0.044	−0.032, −0.770, −0.1545, 0.964, 0.032	0.017, −0.387, −0.064, 0.472, 0.015
6	−0.783, −2.326, −1.77, 5.309, 0.337	−0.028, −0.393, −0.202, 0.988, 0.027	−0.016, −0.901, −0.043, 0.969, 0.023
7	−0.062, −0.483, −0.128, 1.147, 0.039	−0.123, −1.233, −0.315, 1.657, 0.074
8	0.176, −2.467, 0.321, 0.798, 0.111	0.011, −1.164, 0.003, 0.902, 0.001	−0.023, −0.406, −0.123, 0.613, 0.035
9
10	−0.017, −0.990, −0.056, 1.095, 0.016	−0.055, −1.310, −0.188, 1.357, 0.044	−0.016, −0.678, −0.164, 0.963, 0.028
11	−0.039, −2.334, −0.408, 1.666, 0.0467	0.007, −0.968, −0.113, 1.117, 0.005	0.034, −0.778, −0.036, 0.825, −0.009
12	−0.008, −0.926, −0.058, 1.081, −0.002	−0.006, −0.901, −0.071, 0.925, 0.019	−0.018, −0.589, −0.126, 0.785, 0.031
13	−0.069, −1.395, −0.284, 1.565, 0.043	0.000, −1.036, −0.090, 0.948, 0.022	−0.048, −0.526, −0.195, 0.863, 0.050
14	−0.593, −2.030, −0.592, 2.400, 0.309	−0.073, −0.578, −0.237, 1.081, 0.054	−0.012, −0.575, −0.097, 0.707, 0.029
15
16	−0.023, −0.928, −0.100, 1.195, 0.021	−0.045, −1.246, −0.146, 1.301, 0.038

and

Table 9. Surface Reflectance Coefficients for 0.67 µm for DJF.

Land Cover	0.10 < NDVI’ < 0.20	0.20 ≤ NDVI’ < 0.50	0.5 ≤ NDVI’
1	−0.114, −0.704, −0.167, 1.165, 0.083	−0.091, −1.735, −0.198, 1.412, 0.073	−0.086, −0.486, −0.338, 1.188, 0.057
2			−0.387, −0.532, −1.176, 3.248, 0.180
3		0.080, −2.328, 0.056, 1.233, −0.04
4		−0.074, −1.169, −0.192, 1.238, 0.049	−0.377, −0.622, −1.24, 3.304, 0.188
5	0.118, −1.00, 0.138, 0.779, −0.064	−0.082, −1.513, −0.140, 1.386, 0.050	−0.094, −0.513, −0.264, 0.989, 0.065
6
7	0.084, 0.511, −0.483, 1.49, 0.070	−0.196, −1.523, −0.562, 2.072, 0.125	−0.071, −0.522, −0.234, 0.911, 0.060
8	−0.388, −1.127, −0.479, 1.937, 0.203	−0.018, −2.23, 0.011, 1.037, 0.047	−0.070, −0.390, −0.296, 1.008, 0.050
9
10	−0.011, −0.703, −0.058, 1.071, 0.014	−0.057, −1.501, −0.1256, 1.337, 0.036	−0.033, −0.518, −0.157, 0.780, 0.037
11
12	−0.027, −1.456, −0.111, 1.323, 0.017	−0.036, −1.222, −0.193, 1.254, 0.032	−0.090, −0.446, −0.321, 1.065, 0.061
13
14	0.056, −1.074, 0.348, 0.509, −0.034	−0.021, −1.902, −0.075, 1.311, 0.026	0.186, −0.376, 0.148, 0.005, −0.066
15	−0.078, −3.385, −0.034, 1.416, 0.048	0.043, −1.000, −0.114, 1.096, −0.013	−0.009, −0.465, −0.103, 0.597, 0.027
16	−0.070, −0.736, −0.139, 1.291, .0326	−0.041, −1.349, −0.149, 1.311, 0.030	0.018, −0.544, −0.079, 0.655, 0.009

for 0.47 and 0.67 μm GF4 bands, respectively. Corresponding surface types of Land cover values in Table 2, Table 3, Table 4, Table 5, Table 6, Table 7, Table 8 and Table 9 are included in

Table A1. Land Cover types corresponding to the integer in Table 2, Table 3, Table 4, Table 5, Table 6, Table 7, Table 8 and Table 9. Defined according to the International Geosphere-Biosphere Programme (IGBP).

No	Land Cover	No	Land Cover
1	Evergreen Needleleaf Forest	9	Savannas
2	Evergreen Broadleaf Forest	10	Grasslands
3	Deciduous Needleleaf Forest	11	Permanent Wetlands
4	Deciduous Broadleaf Forest	12	Croplands
5	Mixed Forests	13	Urban Areas
6	Closed Shrublands	14	Cropland - Natural Vegetation Mosaic
7	Open Shrublands	15	Snow and Ice
8	Woody Savannas	16	Barren or Sparsely Vegetated

(see Appendix A).

The surface reflectances for selected wavelengths are calculated for each valid observation and subsequently fitted as a function of NDVI’ and scattering angles. This manner allows to capturing the potential seasonal/annual variability of the surface reflectance, enables one to scale this quantity down to the spatial resolution of GF4 and uses it in the retrieval of AOD at GF4 pixel level [26]. At the meantime, the regression has been performed after filtering surface reflectance ‘outliers’ which can be caused by the strong inhomogeneity inside the 5 km land cover resolution. However, the spatial resolution problem still brings larger impact on surface types like urban, where surface variability is large.

Figure 4 shows comparisons of the estimated surface reflectance at 0.67 (a–d) and 0.47 µm (e–h) with the surface reflectance from MOD09A1 data for four seasons: spring (March-April-May, MAM), summer (June-July-August, JJA), autumn (September-October-November, SON), and winter (December-January-February, DJF) in 2009. The color of the symbol indicates the AOD value for each pixel, the AOD is obtained from MODIS C6 DT-DB combined aerosol product [5]. The comparisons show that the estimations for both 0.47 and 0.67 µm have good correlations with the MOD09 surface reflectance for all four seasons, except for high AOD cases. The regression pattern for 0.67 µm is slightly better than 0.47 µm as presented in Figure 4, which may be due to the potential higher chance of aerosol contamination in the atmospheric correction procedure of MODIS surface product for shorter wavelengths [43]. The performance of the atmospheric correction algorithm degrades with increase of AOD values and the atmospheric correction is also less accurate for shorter wavelengths [43]. The regions in Eastern Asia are often covered by high aerosol loading, especially in Eastern China [44]. However, surface reflectance at 0.67 µm for high AOD cases show higher systematical bias compared to 0.47 µm, as the yellow/red dots in Figure 4b and 4c are slightly away from the 1:1 line. This is likely due to the fact that the slopes of surface reflectance from 0.86 to 0.67 μm are more sensitive to the NDVI’ values than the 0.86 to 0.47 μm slopes, as shown in Figure 3. Thus, bias in NDVI’ caused by the presence of aerosols has a stronger impact on the estimated surface reflectance at 0.67 than at 0.47 μm following the above statistical analysis [6].

In addition, we also calculated season-average surface reflectance at 0.47 and 0.67 µm channels for 2009 based on SAHARA regression coefficients listed in Table 2, Table 3, Table 4, Table 5, Table 6, Table 7, Table 8 and Table 9 and DB coefficients [6], respectively. SAHARA and DB estimated surface reflectances are also compared to the MOD09 surface reflectance products in two visible channels as shown in Figure 5 (0.47 µm) and Figure 6 (0.67 µm). All three surface reflectances show similar distribution. The surface reflectances are relatively lower in the Southern part of China compared to the Northern and Western parts for selected wavelengths. At the meantime, the difference between SAHARA surface reflectance and MOD09 is smaller compared to the difference between DB reflectance and MOD09. For example, the overall differences for all seasons are 1.9%, 7.3%, 2.1%, 5.8% for SAHARA and MOD09 for 0.47 µm and 1.2%, −24%, −8%, −6.8% for DB and MOD09, respectively. As to 0.67 µm, the overall difference for all seasons are 8.9%, 5.3%, 8.3%, 6.6% for SAHARA and MOD09 and −15.1%, −16.4%, −22.1%, −18.5% for DB and MOD09. We can also find that the estimation based on DB coefficient tends to underestimate the surface reflectance in most of the seasons and land types causing slight overestimation of AOD over China [44]. The SAHARA coefficients estimated surface reflectances are generally larger than the MOD09 product.

3.2. Aerosol Type Parameteirzation

The aerosol types were predefined according to the season and geo-location, they are: weakly absorbing, moderately absorbing and strongly absorbing, adapted from the current MODIS C6 DT aerosol retrieval algorithm [5]. Levy et al. [41] performed a cluster analysis of the entire AERONET dataset of sun-sky inversions [45] up to 2010. Following this way, an “expected aerosol type” map for all locations at 1° × 1° spatial resolution for each season has been created and used for aerosol retrieval. Figure 7 presents the three absorption aerosol types for the MODIS DT aerosol retrieval algorithm based on the AERONET information [41]. The asymmetry factor (g) and Single Scattering Albedo (SSA) were further parameterized to different aerosol types. After numerical calculation, it is found that the g, SSA and the Angstrom coefficient (Alpha) can be parameterized by the AOD at 0.55 µm using polynomial form. g can be parameterized as follows:

$$g=a_0+a_1{\tau }_{550nm}+a_2{{\tau }^2}_{550nm}$$

(8)

The SSA and Alpha are also parameterized similar as the parameterization for g. The coefficients (a₀, a₁, and a₂) are listed in the

Table 10. The polynomial coefficients for the Single Scattering Albedo (SSA).

SSA	a0	a1	a2
Weakly absorbing	0.921	0.049	−0.018
Moderately absorbing	0.899	0.048	−0.012
Strongly absorbing	0.831	0.044	−0.018

,

Table 11. The polynomial coefficients for the Angstrom coefficient.

Alpha	a0	a1	a2
Weakly absorbing	1.926	−0.217	−0.162
Moderately absorbing	1.865	0.052	−0.268
Strongly absorbing	2.028	0.008	−0.096

and

Table 12. The polynomial coefficients for the asymmetry factor.

g	a0	a1	a2
Weakly absorbing	0.607	0.081	−0.014
Moderately absorbing	0.59	0.053	−0.003
Strongly absorbing	0.548	−0.003	0.024

for SSA, Alpha, and g, respectively.

In order to evaluate the accuracy of the aerosol type parameterization, a comparison between proposed parameterization and AERONET observations was performed. AERONET observations corresponding to time span of GF4 observations for selected regions have been used. Firstly, AODs at 0.55 µm from AERONET level 2.0 were selected from AEROSOL OPTICAL DEPTH (v2) and calculate single scattering albedo and asymmetry factor using Table 10, Table 11 and Table 12, and then the calculated values were compared to AEROSOL INVERSIONS (v2). Figure 8 shows that approximately 85% of estimated asymmetry factor falls into the 10% error envelope (EE) while about 94% single scattering albedo meet the accuracy of 10% EE. About 50% estimated asymmetry factor and 70% estimated SSA fall into 5% EE. Potential overestimation of SSA for strong absorbing aerosol type and a slight underestimation of SSA for moderately/weakly absorbing aerosol types can be found from Figure 8, indicating possible overestimation/underestimation of AOD in SAHARA algorithm depending on aerosol absorption. The comparison with AERONET shows that the aerosol type parameterization provides accurate estimation for aerosol retrieval.

3.3. Aerosol Retrieval Algorithm

The reflectance at the top of the atmosphere (TOA) is described by Equation (9) [12,46]

$$R_{TOA}\left(\lambda ,{\mu }_0,\mu ,\phi \right)=T_g\left[R_{R+A}\left(\lambda ,{\mu }_0,\mu ,\phi \right)+\frac{A\left(\lambda \right)\cdot T\left(\lambda ,{\mu }_0,\mu \right)}{1-A\left(\lambda \right)\cdot s\left(\lambda \right)}\right]$$

(9)

where $\theta =\arccos \mu$ is the satellite zenith angle, ${\theta }_0=\arccos {\mu }_0$ is the solar zenith angle, $\phi$ is the relative azimuth angle, $R_{TOA}\left(\lambda ,{\mu }_0,\mu ,\phi \right)$ is the contribution of the Earth surface and atmosphere to the TOA reflectance, $R_{R+A}\left(\lambda ,{\mu }_0,\mu ,\phi \right)$ is the contribution of the atmosphere reflectance to the TOA reflectance. $T_g$ is the gas absorption transmission, $\mathrm{A}\left(\mathsf{\lambda }\right)$ is the surface spectral albedo, $T\left(\lambda ,{\mu }_0,\mu \right)$ is total transmission and $s\left(\lambda \right)$ is the atmospheric hemispherical albedo and the basic definition of the above parameters can be found below.

Kokhanovsky [47] provided parameterization equations for $s\left(\lambda \right)$ and $T\left(\lambda ,{\mu }_0,\mu \right)$ using the AOD information, which are also used in the SAHARA algorithm. In the SAHARA algorithm, the atmosphere path reflectance is calculated from the Rayleigh reflectance ${\rho }_{Ray}$ and the aerosol reflectance ${\rho }_{aer}$. Equations (10)–(12) describe the calculation of Rayleigh reflectance [48]. $P_{Ray}\left(\Phi \right)$ is the Rayleigh scattering phase function, ${\tau }_{Ray}$ is Rayleigh Optical Depth (ROD) and $\Phi$ is the phase angle between incident direction and scattering direction.

$${\rho }_{Ray}=\frac{{\tau }_{Ray}\cdot P_{Ray}\left(\Phi \right)}{4\mu {\mu }_0}$$

(10)

$${\tau }_{Ray}=0.00864{\lambda }^{-\left(3.916+0.074+\frac{0.05}{\lambda }\right)}$$

(11)

$$P_{Ray}\left(\Phi \right)=\frac{3}{4}\left(1+{\cos }^2\Phi \right)$$

(12)

The aerosol reflectance can be parameterized to the aerosol Single Scattering Albedo (SSA, ϖ), aerosol phase function ($P_{aer}\left(\Phi \right)$) and AOD ($\mathsf{\tau }$). The single scattering approximation can be used for relatively small AOD, which has been used in the Aerosol Retrieval Technique (ART) [10] and SARA algorithm [7]. Equation (13) describes the single scattering approximation of the aerosol reflectance. The Henyey-Greenstein (H-G) aerosol phase function is used [49], as Equation (14) describes, where g is the asymmetry factor.

$${\rho }_{aer}=\frac{\varpi \cdot \tau \cdot P_{aer}\left(\Phi \right)}{4\mu {\mu }_0}$$

(13)

$$P_{aer}\left(\Phi \right)=\frac{1-g^2}{{(1+g^2+2g\cos \Phi )}^{\frac{3}{2}}}$$

(14)

Then the optimal AODs are obtained by minimize the difference of TOA reflectance at 0.47 and 0.67 µm between observed by GF4 and calculated by Equation (9) using Levenberg–Marquardt iteration algorithm [34]. Figure 9 show the flowchart of SAHARA algorithm.

4. Results and Discussion

In order to show the accuracy of the SAHARA algorithm, the complicated scenes over Eastern and Middle China are selected for testing. Please note that GF4 starts providing stable observations with relative good calibration accuracy in May 2016, the calibration radiance accuracy of GF-4 PMS images are better than 3% [50]. Figure 10 and Figure 11 show the comparisons of AOD retrieved by SAHARA using GF4 data and MODIS C6 DT-DB combined AOD over different surface types. Considering the urgent requirement of high spatial resolution AOD for the investigation of aerosol properties in a local scale, the figures contain the 10 km AOD results for both MODIS and GF4 and 50 meter GF4 retrieved AODs to illustrate the capability of the usage of SAHARA algorithm to very high spatial resolution aerosol retrievals.

Figure 10 presents the 10 km (a and b) and 50 meter (c) resolution AOD at 0.55 µm over Shandong Peninsula and Jiangsu province on 25 June 2016, from MODIS (a) and SAHARA retrieval (b and c), respectively. The spatial distribution of AOD shows large variability for this scenario with minimal and maximal AOD about 0.15 and 1.5. Both MODIS and GF4 retrievals catch the large variability. MODIS and GF4 results are slightly different for the high AOD patterns, mainly due to the different overpass time. The overpass time for GF4 and MODIS/Aqua, used for comparison are 07:00 UTC and 05:10 UTC. According to the wind data from the dataset (V3.0) of daily values of climate data from Chinese surface stations for global exchange, the daily-average wind-speed over Shandong Peninsula is about 5 m/s, and the direction of the maximum wind speed is southwest on 25 June 2016. The southwest wind blows the polluted air mass to the northeast, causing higher aerosol loading over the north of Shandong peninsula in SAHARA retrieval compared to MODIS retrieval. According to the two 10 km retrievals, we can see that the AODs for the upper part of the image (Shandong Peninsula) are much smaller compared to the lower part (Jiangsu Province). The main reason for this is that there are more industries in Jiangsu province while the main aerosol source because the straw combustion season (starting from late August) over Shandong has not yet started [51]. We can clearly see a large difference near the cloud edge (lower part of the image), indicating potential cloud contamination for GF4 aerosol retrieval. When comparing the 50 meter AOD to the 10 km resolution AOD, a good local transportation pattern can been seen from the 50 meter result while this is not the case for 10 km results. Another evidence of possible cloud contamination is that there are some red outliners in GF4 10 km while this is not the case for 50 meter AOD because potential cloud contributes to the averaged 10 km GF data. Another interesting thing is that GF4 10 km retrievals lose information near the coastline while the 50 meter AOD has more coverage close to the water. A similar pattern can be found for Figure 11, illustrating the AOD over Inner Mongolia, which is mainly covered by grasslands. The shift of the high AOD pattern can be explained by the difference of overpass time, the GF-4 image is obtained at 03:00 UTC, while the MODIS result is from 04:10 UTC. The meteorological condition can be another possible reason for the quick change of the aerosol spatial distribution between GF4 and MODIS. During this day, west wind is prevailing over Inner Mongolia, the maximum speed reached 10 m/s and the average speed is about 5–6 m/s.

Figure 10 and Figure 11 also show the dynamical presentation of the parameterization method proposed in the paper for different months, which is linked to the seasonal/geo-location dependent aerosol types. For the current retrieval, a fixed AOD a-priori value (0.2), which approaches the AOD global mean value, is used [4]. According to Figure 10 and Figure 11, the SHARA AOD agrees quite well with MODIS C6 aerosol product. However, the SAHARA AODs are lower compared to MODIS product in the high aerosol loading regions. This can be explained from four aspects. Firstly, the MODIS aerosol product over China is validated to be overestimated [44,52]. Secondly, the parameterization method proposed in this paper works pretty well for AODs range from 0.1 to 1, which may produce some uncertainty for large AODs [53,54]. Thirdly, the underestimation of SSA for aerosol type parameterization may lead to the underestimation of AOD. Fourthly, the SAHARA surface parameterization may slightly overestimate the surface reflectance cause underestimation of AOD. We also notice that there are higher chances to retrieve negative values due to the calibration uncertainty in GF4 compared to MODIS observations. The negative values in MODIS aerosol product [41] indicate the potential calibration error and the uncertainty of calibration error can be enhanced due to the simplification in SAHARA algorithm, which also increases with the increase of the aerosol absorption characteristics.

The case studies show that the SAHARA algorithm is able to capture the aerosol information, and can obtain high resolution AOD. We perform an additional comparison with MODIS aerosol product and validate using AERONET observations. We quantitatively compare SAHARA-derived 10 km AOD and MODIS C6 combined AOD in Figure 12b for the scenario shown in Figure 10. SAHARA AOD shows good agreement with MODIS AOD, with a correlation coefficient of R² = 0.72. All GF4 observations starting from May 2016 are used in the validation with AERONET observations. Due to limited GF4 data and ground-based observations (all available and released GF4 data from 16 May to 30 August have been used), we can only obtain several station data after the collocation. Figure 12a shows the scattering plot of AERONET AOD and SAHARA 50 meter AOD. Good agreement between SAHARA AOD and AERONET AOD, with correlations coefficient R² = 0.86. Besides, the Root-Mean-Square Error (RMSE) is about 0.13 and the Mean Absolute Error (MAE) is 0.12. About 78% of the SAHARA retrieved AOD fall within the EE (Expected Error, $0.20\ast AO{D}_{AERONET}\pm 0.05$) envelope ($(AO{D}_{AERONET}-\left|EE\right|)<AO{D}_{SAHARA}<(AO{D}_{AERONET}+\left|EE\right|)$) [5], which indicates the retrieved AOD of good quality.

5. Conclusions

A new aerosol retrieval algorithm, SAHARA is presented, which enables the retrieval of the AOD, which will be further used for the atmospheric correction for Chinese GaoFen4 instrument. Without using LUT, the SAHARA algorithm uses the parameterized radiative transfer equation for different aerosol types. They are weakly absorbing, moderately absorbing and strongly absorbing, and are adapted from the current MODIS DT C6 aerosol retrieval algorithm [5]. The aerosol types are predefined according to season and geo-location. The asymmetry factor, SSA and the Angstrom coefficient are parameterized by the AOD at 0.55 µm using polynomial form for each aerosol type, which is proved to be accurate enough for the AOD range from 0.1 to 1.0. SAHARA algorithm has a surface parameterization following a similar idea to the MODIS DB algorithm. The surface reflectances at GF4 visible channels have been parameterized by the Rayleigh-corrected TOA reflectance at 0.86 µm, vegetation amount determined by the NDVI’ as well as the scattering angles. Seasonal based coefficients for different surface types are obtained by fitting to the MODIS surface reflectance product. Comparison between SAHARA, DB surface parameterization and MODIS surface reflectance products shows good agreement. SAHARA parameterization shows slight overestimation of surface reflectance while DB estimations underestimate compared to MODIS surface reflectance product. Both SAHARA and DB surface parameterizations show high accuracy for aerosol retrieval.

Two complicated scenes over China are used to test the SAHARA algorithm, the retrieval results from SAHARA agree quite well with the MODIS aerosol product for 10 km spatial resolution. SAHARA is applied for the retrieval of GF4 data based on the original pixel resolution (50 meter) and the retrieval results of 50 meter show reasonable patterns for aerosol distribution. The transportation features of high aerosol loading are well-caught. The validation using AERONET observations shows good performance with correlation coefficient of R² = 0.86 and the slope of 1.05 shows acceptable accuracy of the parameterization of aerosol types and relatively small bias (0.01), showing good accuracy of surface parameterization.

SAHARA has the potential ability to retrieve AOD over both bright and dark surfaces. The retrieved AOD can be used for atmospheric correction of GF4 data which have a spatial resolution of 50 meter. This algorithm can also be used in other satellite data which have similar bands with GF4 instruments, such as Landsat and Sentinel observations [55]. The retrieved AOD shows great potential in monitoring the local air pollution with high spatial resolution to capture local pollution sources. The retrieved AOD can also be used for data assimilation with other satellite data (MODIS, Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) aerosol product, etc.), ground-based measurements (AERONET, particulate matter (PM) observations, and meteorological observations) and model data (GOCART aerosol data) to further improve the air quality analyses and forecasts [56,57].

There are several aspects to be considered regarding to the SAHARA aerosol retrieval algorithm in the future. The first aspect is the cloud screening. The MODIS DT cloud screening algorithm without cirrus cloud test (1.38 µm test) has been used with the current retrieval together with a visual-check according to the RGB images. An individual team is still working on the development of cloud detection of GF4, which will be used in the next retrieval version of SAHARA. The second aspect is linked to the surface parameterization. Currently only one year data were used for the statistical analysis and there are no coefficients for certain surface types of some seasons, which will be solved by extending the regression dataset. The possible uncertainty of surface reflectance estimation due to the BRDF effect, especially for high spatial resolution AOD retrieval, needs to be considered. The Ross-Li BRDF model will be implemented in future work [58,59]. The third aspect is linked to the aerosol type parameterization. According to previous evaluation, the MODIS DT algorithm provides overestimated AOD partly due to the aerosol type, which is parameterized by limited AERONET observation. Thus, an updated aerosol type parameterization will be necessary based on both AEROENT and CARSNET and other in-situ measurements.