A Dual-Channel Aerosol Optical Depth Retrieval Algorithm Incorporating the BRDF Effect from AVHRR over Eastern Asia

: A NOAA/AVHRR dual-channel method over land is proposed to simultaneously retrieve aerosol optical depth (AOD) at 0.55 µ m, and surface reﬂectance at 0.63 and 0.85 µ m. Compared with previous well-established one-channel retrieval algorithms, this algorithm takes advantage of the surface reﬂectance ratio between 0.63 and 0.85 µ m in an attempt to account for the effect induced by the surface bidirectional reﬂectance distribution function (BRDF). This effect cannot be negligible due to the orbit drift and phasing running of NOAA satellites, both of which potentially cause large solar angular variation. Meanwhile, the observation posture change of AVHRR would cause large sensor angular variation in time series measurements. The used surface reﬂectance ratio based on dual channels at 0.63 and 0.85 µ m is found more reasonable to be assumed as unchanged during a certain period of time, compared to the traditional ratio when addressing the BRDF issue. AOD retrievals have been carried out over Eastern Asia. Validation against aerosol robotic network (AERONET) measurements shows that up to 83% of AOD validation collocations are within error lines ( ± 0.15 ± 0.15 τ , τ is AOD) with an R of 0.88 and an root mean square error (RMSE) of 0.15. The dual-channel algorithm taking into account the surface BRDF effect is proved outperforming the conventional 0.63 µ m-channel method. It indicates that our algorithm has the potential to be applied to the retrieval of long series of AOD, especially to the AOD retrieval of the sensors which lack a shortwave infrared channel required in the MODerate resolution Imaging Spectroradiometer (MODIS) dark target AOD algorithm. BRDF effect in the determination of surface reﬂectance from time series of AVHRR measurements. Our forward sensitivity study found the need for dual channels at 0.63 µ m and 0.85 µ m to retrieve AOD instead of the single-channel algorithm using 0.63 µ m. We used the surface reﬂectance ratio k (cid:48) between 0.63 µ m and 0.85 µ m to address the surface BRDF effect. The ratio was demonstrated from both the forward-sensitivity and theoretical point of view and the validity was examined. It was found that the AOD inversion using the k (cid:48) ratio is more robust to the input error than the conventional 0.63 µ m-channel algorithm. Moreover, it is more reasonable to assume that the k (cid:48) ratio remains unchanged during a certain period of time compared to the traditional ratio k when addressing the BRDF issue. Then, the AVHRR dual-channel algorithm was proposed to determine AOD at 0.55 µ m and surface reﬂectance at 0.63 and 0.85 µ m. This novel approach was applied to Eastern Asia. The surface reﬂectance retrieval was compared with the atmospherically corrected surface reﬂectance value and MODIS surface reﬂectance product. The surface reﬂectance at 0.85 µ m has better statistic results than that at 0.63 µ m, consistent with the uncertainty conclusion. The AOD retrieval was compared with AERONET AOD observation, yielding the linear regression equation of y = 0.074 + 0.89 x with an R of 0.88, an RMSE error of 0.15. The fraction of the AOD validation collocations within the error lines ( ± 0.15 ± 0.15 τ ) was 83%. The favorable correlation provides conﬁdence in the utility of our AVHRR AOD algorithm over land. Time series comparison between AVHRR AOD retrievals and AERONET shows that the AVHRR AOD retrievals can capture the AOD changes over time. dual-channel AOD


Introduction
Aerosols play an important role in the regional and global climate through direct interactions with radiation [1][2][3][4] and indirect interactions with clouds [5][6][7][8][9][10][11]. It has been well recognized that satellite retrieved aerosol optical depth (AOD) data is one of the useful data sources to constrain the uncertainties of climate effect induced by intensive anthropogenic activities (e.g., [12]). In particular, the long-term accurate AOD measurement, from either ground-based or satellite-based instruments, can well reflect the trend of local and global air pollution [7,13,14].
Ground-based AOD observations from the aerosol robotic network (AERONET) can provide AOD measurements with high accuracy [15]. However, AERONET is not available before 1997 in eastern Asia. Moreover, it only reflects local aerosol properties near the ground-based stations. In comparison, the Advanced Very High-Resolution Radiometer in the imaging time of the same pixel on adjacent days reaching ±51 min, leading to the solar angular variation unignorable [35,36]. Besides, the large AVHRR observing posture variation causes the large variation in view geometry. All the above factors lead to an unignorable BRDF effect. Only a few studies have accounted for the BRDF effect [19,20,22]. Further, the surface reflectance at 0.63 µm alone is used to retrieve AOD. However, our forward sensitivity study reveals that in some cases only using the 0.63 µm channel to retrieve AOD may be inadequate and cause large error. Therefore, a dual-channel retrieval algorithm is required. Thus, we go further to take advantage of the surface reflectance ratio between 0.63 and 0.85 µm to estimate the clear day (the least TOA reflectance ratio) and to retrieve AOD instead of the surface reflectance at 0.63 µm. This paper is organized in the following manner. Section 2 introduces the data sets and the preprocess method. The whole algorithm is depicted in detail in Section 3. The need of a dual-channel retrieval algorithm is described firstly. The reason for choosing the surface reflectance ratio to address the BRDF effect is then demonstrated from the theoretical point of view, including the uncertainty study as well. The validation results are presented in Section 4. Finally, the key conclusions are summarized in Section 5.

Datasets and Their Preprocessing
We mainly used AVHRR observations as the data source. Besides, AERONET measurements were used to analyze the dual-channel retrieval algorithm and validate the AOD retrieval result. The two datasets are described below.

AVHRR
The Local Area Coverage (LAC) data of AVHRR/3 was used in this study. AVHRR/3 observes the Earth in the following six spectral bands: 0.58-0.68 µm (channel 1), 0.725-1.0 µm (channel 2), 1.58-1.64 µm (channel 3A), 3.55-3.93 µm (channel 3B), 10.3-11.3 µm (channel 4), and 11.5-12.5 µm (channel 5). The nominal resolution of the LAC data at the nadir is 1.1 km. NOAA uses two satellites, namely a morning satellite and an afternoon satellite, to ensure that the sunlight side of the Earth is observed at least twice every day. In this study, Eastern Asia was selected as the study area ( Figure 1). The geographic range covers from 18 • N to 54 • N and 110 • E to 136 • E. AVHRR data from July 2008 to June 2010 were collected. During this period, three satellites were involved. NOAA-18 (NN) and NOAA-19 (NP) are the afternoon satellites, while NOAA-17 (NM) is the morning satellite.
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 21 be even larger for the afternoon satellites. Moreover, the phasing running of the NOAA satellites can result in the imaging time of the same pixel on adjacent days reaching ±51 min, leading to the solar angular variation unignorable [35,36]. Besides, the large AVHRR observing posture variation causes the large variation in view geometry. All the above factors lead to an unignorable BRDF effect. Only a few studies have accounted for the BRDF effect [19,20,22]. Further, the surface reflectance at 0.63 μm alone is used to retrieve AOD. However, our forward sensitivity study reveals that in some cases only using the 0.63 μm channel to retrieve AOD may be inadequate and cause large error. Therefore, a dual-channel retrieval algorithm is required. Thus, we go further to take advantage of the surface reflectance ratio between 0.63 and 0.85 μm to estimate the clear day (the least TOA reflectance ratio) and to retrieve AOD instead of the surface reflectance at 0.63 μm. This paper is organized in the following manner. Section 2 introduces the data sets and the preprocess method. The whole algorithm is depicted in detail in Section 3. The need of a dual-channel retrieval algorithm is described firstly. The reason for choosing the surface reflectance ratio to address the BRDF effect is then demonstrated from the theoretical point of view, including the uncertainty study as well. The validation results are presented in Section 4. Finally, the key conclusions are summarized in Section 5.

Datasets and Their Preprocessing
We mainly used AVHRR observations as the data source. Besides, AERONET measurements were used to analyze the dual-channel retrieval algorithm and validate the AOD retrieval result. The two datasets are described below.

AVHRR
The Local Area Coverage (LAC) data of AVHRR/3 was used in this study. AVHRR/3 observes the Earth in the following six spectral bands: 0.58-0.68 μm (channel 1), 0.725-1.0 μm (channel 2), 1.58-1.64 μm (channel 3A), 3.55-3.93 μm (channel 3B), 10.3-11.3 μm (channel 4), and 11.5-12.5 μm (channel 5). The nominal resolution of the LAC data at the nadir is 1.1 km. NOAA uses two satellites, namely a morning satellite and an afternoon satellite, to ensure that the sunlight side of the Earth is observed at least twice every day. In this study, Eastern Asia was selected as the study area ( Figure 1). The geographic range covers from 18°N to 54°N and 110°E to 136°E. AVHRR data from July 2008 to June 2010 were collected. During this period, three satellites were involved. NOAA-18 (NN) and NOAA-19 (NP) are the afternoon satellites, while NOAA-17 (NM) is the morning satellite. The AVHRR data preprocess and gas absorption correction method is described as follows. Firstly, the AVHRR TOA radiance was calibrated into reflectance using the postlaunched calibration coefficients from NOAA/National Environmental Satellite The AVHRR data preprocess and gas absorption correction method is described as follows. Firstly, the AVHRR TOA radiance was calibrated into reflectance using the post-launched calibration coefficients from NOAA/National Environmental Satellite Information Data and Information Service (NOAA/NESDIS) [37,38]. Secondly, the clouds were masked using the Clouds from AVHRR Extended System (CLAVR-x) results [39]. We also filtered the oceanic pixels using MODIS land water mask product (MOD44W). Followingly, because snow-or ice-covered surfaces reflect too much light, snow and ice pixels were eliminated, especially in the winter and spring seasons of Eastern Asia. The snow and ice product from global map datasets at 4 × 4 km spatial resolution (ftp://www.orbit.nesdis.noaa.gov/pub/smcd/emb/snow/global_mult_snow_ice/) were used [40]. Then, pixels of which satellite zenith angles larger than 60 • or solar zenith angles larger than 63 • were discarded because plane parallel approximation was used in the study.
The AVHRR TOA reflectance at 0.63 and 0.85 µm should be corrected to remove the effects of water vapor and ozone because 0.85 µm channel is heavily influenced by a strong water vapor absorption band. The MODIS gas correction method was employed (http://modis.gsfc.nasa.gov/data/atbd/atmos_atbd.php). Products from the Ozone Monitoring Instrument (OMI) on board the NASA Aura satellite were used as ozone input. We recalculated the water vapor absorption coefficients which were adjusted for AVHRR spectral responses using the Second Simulation of a Satellite Signal in the Solar Spectrum-Vector (6SV) radiative transfer code [41]. The formula adopted for water vapor transmission is where M is the air mass, U H 2 O is the total precipitable water vapor and we used MOD/MYD07 products as input. Coefficients a, b, and c were calculated, and the values of a, b, and c are listed in Table 1. The values of a, b, and c in 0.63 µm change little from those in MODIS 0.66 µm band, while a, b, and c values in 0.85 µm are significantly different from those in MODIS 0.86 µm band. It is revealed that different sensor spectral responses ( Figure 2) over the water vapor absorption band can cause different water vapor influences. Therefore, the water vapor correction of AVHRR 0.85 µm channel is necessary. In the absence of MOD/MYD07 water vapor data, the global average optical depth values were used and the Beer-Bouguer-Lambert law was adopted (http://modis.gsfc.nasa.gov/data/atbd/atmos_atbd.php). Information Data and Information Service (NOAA/NESDIS) [37,38]. Secondly, the clouds were masked using the Clouds from AVHRR Extended System (CLAVR-x) results [39]. We also filtered the oceanic pixels using MODIS land water mask product (MOD44W). Followingly, because snow-or ice-covered surfaces reflect too much light, snow and ice pixels were eliminated, especially in the winter and spring seasons of Eastern Asia. The snow and ice product from global map datasets at 4 × 4 km spatial resolution (ftp://www.orbit.nesdis.noaa.gov/pub/smcd/emb/snow/global_mult_snow_ice/) were used [40]. Then, pixels of which satellite zenith angles larger than 60° or solar zenith angles larger than 63° were discarded because plane parallel approximation was used in the study.
The AVHRR TOA reflectance at 0.63 and 0.85 μm should be corrected to remove the effects of water vapor and ozone because 0.85 μm channel is heavily influenced by a strong water vapor absorption band. The MODIS gas correction method was employed (http://modis.gsfc.nasa.gov/data/atbd/atmos_atbd.php). Products from the Ozone Monitoring Instrument (OMI) on board the NASA Aura satellite were used as ozone input. We recalculated the water vapor absorption coefficients which were adjusted for AVHRR spectral responses using the Second Simulation of a Satellite Signal in the Solar Spectrum-Vector (6SV) radiative transfer code [41]. The formula adopted for water vapor transmission is where is the air mass, is the total precipitable water vapor and we used MOD/MYD07 products as input. Coefficients , , and were calculated, and the values of , , and are listed in Table 1. The values of , , and in 0.63 μm change little from those in MODIS 0.66 μm band, while , , and values in 0.85 μm are significantly different from those in MODIS 0.86 μm band. It is revealed that different sensor spectral responses ( Figure 2) over the water vapor absorption band can cause different water vapor influences. Therefore, the water vapor correction of AVHRR 0.85 μm channel is necessary. In the absence of MOD/MYD07 water vapor data, the global average optical depth values were used and the Beer-Bouguer-Lambert law was adopted (http://modis.gsfc.nasa.gov/data/atbd/atmos_atbd.php).   After the gas absorption correction, the remaining AVHRR pixels were further processed to remove undetected clouds or cloud shadow pixels and resized to a resolution Remote Sens. 2021, 13, 365 5 of 21 of 11 km, because cloud shadow may contaminate the surface reflectance and mistaking clouds for aerosol could lead to abnormal high AOD retrievals. The brightest 50% and darkest 20% of AVHRR pixels were removed when the TOA reflectance values of channel 1 were between 0 and 0.25, which is similar as what is done in the MODIS pixel pre-selection method.

AERONET Measurements
AERONET is a global network that provides long-term and continuous AOD data with relatively high accuracy and precision [15,42] for validating satellite-based AOD retrievals [43,44]. Overall, six AERONET sites were chosen to validate the following AVHRR retrievals: Beijing, Ussuriysk, EPA-NCU, Gwangju_GIST, XiangHe, and Xinglong ( Figure 1). Their latitudes, longitudes and elevations are displayed in Table 2. The last column of Table 2 also gives the land surface coverage types, which are roughly outlined in the view of the site information from the AERONET website. In this study, level 2.0 cloud-screened and quality-assured AOD data of version 3 algorithm [15,42] were collected. To compare with the AVHRR-retrieved AOD values at 0.55 µm, ground-based AOD values at 0.55 µm were logarithmically calculated using instantaneous AERONET AOD estimates at 0.675 µm and the Ångström exponent (0.44-0.87 µm). The AERONET sun photometers also sample the almucantar sky radiance every 60 min. The inversion algorithm retrieves the aerosol optical properties, namely size distribution and refractive index, from diffuse sky observations. Level 1.5 size distribution and complex refractive indices products were used for high availability.

Retrieval Methodology
The flowchart of the dual-channel retrieval algorithm is shown in Figure 3. A brief description of the workflow is given as below.
(a) AVHRR data preprocess: The detailed description of calibration, the screening, and the gas absorption correction has been given in Section 2.1. Then time series of AVHRR TOA reflectance at 0.63 and 0.85 µm were prepared. (b) Dual-channel retrieval algorithm: Firstly, the clear day in 31 days was searched for each pixel according to the least TOA reflectance ratio. Secondly, the background AOD and Rayleigh scattering of the clear day was removed. Next, assuming that the surface bidirectional reflectance properties do not change during the 31 days, the surface reflectance ratio of the retrieval day was obtained. Finally, the AOD values at 0.55µm and the surface reflectance values at 0.63 and 0.85 µm were retrieved using the surface reflectance ratio of the retrieval day on the basis that the surface reflectance ratio would be monotonically decreasing as the AOD value increases. Detailed description of the algorithm is shown in Section 3.4.

(c) Comparison and validation:
The retrieved AOD and surface reflectance results are evaluated in terms of statistical comparisons and spatial analysis of one-day case (shown in Section 4). This part begins with the forward sensitivity analysis and the need for a dual-channel algorithm is derived. After that, the reason for choosing the surface reflectance ratio to address the BRDF effect is demonstrated theoretically in Section 3.2. Next, the whole dualchannel retrieval algorithm is stated in detail in Sections 3.3 and 3.4. Then the uncertainty analysis was studied in Section 3.5.

Forward Sensitivity Study and the Need for a Dual-Channel AOD Algorithm
The TOA reflectance values at 0.63 and 0.85 µm were simulated as a function of AOD and the surface reflectance values at 0.63 and 0.85 µm using the 6SV software [41]. The two sets of the surface reflectance values at 0.85 µm (0.15 and 0.2) were used, each one of which was assumed with two different surface reflectance ratios between 0.63 and 0.85 µm (k = 0.1177, 0.845). The solar zenith angle is 22.14 • , the satellite zenith angle is 48 • , and the relative azimuth angle −156.41 • . Figure 4 shows that under some condition, only using 0.63 µm channel to retrieve AOD may cause larger error than using both channels at 0.63 and 0.85 µm. For example, when k = 0.1177, two TOA reflectance lines of 0.63 µm are too close to each other so that it is difficult to accurately extract AOD value from the TOA reflectance values, especially when AOD value is high. To depict the surface reflectance ratio k in such condition, we simulated the TOA reflectance ratio of 0.63 to 0.85 µm as the function of AOD using 6SV software [41] (the relationship is shown in Figure 5). Overall, four groups of surface reflectance values at 0.63 and 0.85 µm were used at the same values of Figure 4. The sun-view angles were also taken the same values. It can be shown that using k instead of the surface reflectance at 0.63 µm can help differentiate two indistinguishable lines in Figure 4 when k = 0.1177. The two lines can be distinguished more easily as the AOD value increases. Therefore, it is necessary to use the measurements of two AVHRR red and near-infrared channels in order to retrieve accurate AOD values. It should also be noted that in Figure 4 when surface reflectance value at 0.85 µm is 0.2 the TOA reflectance value decreases with the increasing AOD. This is contrary to the assumption based on which the clear day is selected and should be avoided in the AOD retrieval. the measurements of two AVHRR red and near-infrared channels in order to retrieve accurate AOD values. It should also be noted that in Figure 4 when surface reflectance value at 0.85 μm is 0.2 the TOA reflectance value decreases with the increasing AOD. This is contrary to the assumption based on which the clear day is selected and should be avoided in the AOD retrieval.  In order to further illustrate the different sensitivities of AOD inversion between using the conventional 0.63 μm channel and using the ′ ratio in the case of ′ = 0.1177, the AOD values at 0.55 μm were retrieved with errors added to the surface reflectance at 0.63 μm and to the ′ ratio, respectively. Errors ranging from −0.05 to 0.05 stepped by 0.01 were used. In the forward simulation the surface reflectance at 0.85 μm was set as 0.15 and the ′ value of 0.1177 was used. The sun-view angles were the same as Figure 4 and 5. It is shown that the AOD mean absolute error (MAE) of the conventional algorithm is much larger than that of the dual-channel algorithm ( Figure 6). Therefore, the AOD retrieval algorithm using the ′ ratio is more robust to the input error than the conventional 0.63 μm-channel algorithm. The founding illustrates that under some the measurements of two AVHRR red and near-infrared channels in order to retrieve accurate AOD values. It should also be noted that in Figure 4 when surface reflectance value at 0.85 μm is 0.2 the TOA reflectance value decreases with the increasing AOD. This is contrary to the assumption based on which the clear day is selected and should be avoided in the AOD retrieval.  In order to further illustrate the different sensitivities of AOD inversion between using the conventional 0.63 μm channel and using the ′ ratio in the case of ′ = 0.1177, the AOD values at 0.55 μm were retrieved with errors added to the surface reflectance at 0.63 μm and to the ′ ratio, respectively. Errors ranging from −0.05 to 0.05 stepped by 0.01 were used. In the forward simulation the surface reflectance at 0.85 μm was set as 0.15 and the ′ value of 0.1177 was used. The sun-view angles were the same as Figure 4 and 5. It is shown that the AOD mean absolute error (MAE) of the conventional algorithm is much larger than that of the dual-channel algorithm ( Figure 6). Therefore, the AOD retrieval algorithm using the ′ ratio is more robust to the input error than the conventional 0.63 μm-channel algorithm. The founding illustrates that under some In order to further illustrate the different sensitivities of AOD inversion between using the conventional 0.63 µm channel and using the k ratio in the case of k = 0.1177, the AOD values at 0.55 µm were retrieved with errors added to the surface reflectance at 0.63 µm and to the k ratio, respectively. Errors ranging from −0.05 to 0.05 stepped by 0.01 were used. In the forward simulation the surface reflectance at 0.85 µm was set as 0.15 and the k value of 0.1177 was used. The sun-view angles were the same as Figures 4 and 5. It is shown that the AOD mean absolute error (MAE) of the conventional algorithm is much larger than that of the dual-channel algorithm ( Figure 6). Therefore, the AOD retrieval algorithm using the k ratio is more robust to the input error than the conventional 0.63 µm-channel algorithm. The founding illustrates that under some condition using the 0.63 µm channel alone to retrieve AOD may be inadequate and cause large error. The effect of the k ratio error on the discrimination of the TOA reflectance ratios was also investigated. The sun-view angles and the surface reflectance values at 0.85 µm were taken the same values as Figures 4 and 5. The error of ±0.05 and ±0.1 was added to the k value (depicted in Figures 7 and 8, respectively). It is shown that for k = 0.845 the error of ±0.05 would cause the TOA reflectance ratio indistinguishable, while for k = 0.1177 the condition varies with the AOD value. When the AOD value is less than 1, the error of ±0.05 would cause the TOA reflectance ratio difficult to be distinguished (Figure 7). When the AOD value is larger than 1, the error of the k value causing the TOA reflectance ratio indistinguishable would be 0.1 (Figure 8). condition using the 0.63 μm channel alone to retrieve AOD may be inadequate and cause large error. The effect of the ′ ratio error on the discrimination of the TOA reflectance ratios was also investigated. The sun-view angles and the surface reflectance values at 0.85 μm were taken the same values as Figures 4 and 5. The error of ±0.05 and ±0.1 was added to the ′ value (depicted in Figures 7 and 8, respectively). It is shown that for ′ = 0.845 the error of ±0.05 would cause the TOA reflectance ratio indistinguishable, while for ′ = 0.1177 the condition varies with the AOD value. When the AOD value is less than 1, the error of ±0.05 would cause the TOA reflectance ratio difficult to be distinguished ( Figure  7). When the AOD value is larger than 1, the error of the ′ value causing the TOA reflectance ratio indistinguishable would be 0.1 (Figure 8).

Why Choosing the Surface Reflectance Ratio to Address the BRDF Effect?
Beside the above forward sensitivity study, the reason for choosing the surface reflectance ratio to address the BRDF effect due to the view geometry variation and solar angular change is also demonstrated from the theoretical point of view followingly.
Previous studies have shown that the surface bidirectional reflectance properties can be described by two parts that vary separately with wavelength and geometry [45]. Moreover, the scale of surface geometrical structure considering the shadowing effects is much larger than the wavelength [45]. Therefore, the surface reflectance shape can be assumed as independent of the wavelength [46]. Furthermore, the surface bidirectional reflectance , can be described as the product of the surface reflectance magnitude and the surface reflectance shape [47].
where , is the surface reflectance magnitude varying with the wavelength and time, is the surface reflectance shape mainly depending on the geometry Ω.
can be assumed as varying much more slowly than the surface reflectance magnitude [47].
Then the following ratio was used [48,49] to address the surface reflectance issue in AOD retrieval as follows.
Herein, Equation (3) is transformed to more conveniently relate the surface bidirectional reflectance of the retrieval day to that of the clear day. The transformation is From Equations (2) and (3) Because is assumed as varying much more slowly than the surface reflectance magnitude [47], Equation (5) When and take the red and near-infrared bands, respectively, the quantity of − can be considerable, as for the soybeans, savannah and pasture land (Figure 9).

Why
Choosing the Surface Reflectance Ratio to Address the BRDF Effect?
Beside the above forward sensitivity study, the reason for choosing the surface reflectance ratio to address the BRDF effect due to the view geometry variation and solar angular change is also demonstrated from the theoretical point of view followingly.
Previous studies have shown that the surface bidirectional reflectance properties can be described by two parts that vary separately with wavelength and geometry [45]. Moreover, the scale of surface geometrical structure considering the shadowing effects is much larger than the wavelength [45]. Therefore, the surface reflectance shape can be assumed as independent of the wavelength [46]. Furthermore, the surface bidirectional reflectance A Ω,λ can be described as the product of the surface reflectance magnitude and the surface reflectance shape [47].
where f λ, t is the surface reflectance magnitude varying with the wavelength λ and time, B Ω is the surface reflectance shape mainly depending on the geometry Ω. B Ω can be assumed as varying much more slowly than the surface reflectance magnitude [47]. Then the following ratio was used [48,49] to address the surface reflectance issue in AOD retrieval as follows.
Herein, Equation (3) is transformed to more conveniently relate the surface bidirectional reflectance of the retrieval day to that of the clear day. The transformation is From Equations (2) and (3), we can obtain Because B Ω is assumed as varying much more slowly than the surface reflectance magnitude [47], Equation (5) can be simplified as When λ 1 and λ 2 take the red and near-infrared bands, respectively, the quantity of k λ 1 − k λ 2 can be considerable, as for the soybeans, savannah and pasture land ( Figure 9). However, if we substitute the denominator of the right part of Equation (6) for f λ 2 ,t 1 · f λ 2 ,t 2 , we can obtain Figure 9. The mean absolute errors of k Ω 1 − k Ω 2 and k λ 1 − k λ 2 for different ground cover types. Equation (7) is generally smaller than Equation (6) because f λ 2 ,t 1 is usually much larger than f λ 1 ,t 2 . Therefore, it is more reasonable to assume that the surface reflectance ratio k remains unchanged during a certain period of time compared to k λ (Equation (3)).
Followingly, the Rahman-Pinty-Verstraete (RPV) model [50] was used to verify the validity of Equation (8). The RPV model was developed through taking the hot spot effect into consideration on the basis of earlier empirical bidirectional reflectance models. A total of three input parameters are required in the model. Rahman et al. [50] has listed a set of inversion results of the three parameters for different cover types in the AVHRR red and near-infrared channel against various datasets, including ground-based, laboratory, and airborne measurements ( Table 1 in the RPV paper). The inversion results were used in the analysis. The hard wheat type was excluded for possible printing error because the surface reflectance values of the hard wheat are different from those in the original data source [51]. The validity of Equation (8) was analyzed as follows. Firstly, the geometries of the test were from the same location at which 15 AVHRR images from 14 June 2009 to 28 June 2009 onboard the NOAA-19 satellite were taken. Secondly, the surface reflectance values were calculated from the RPV model using the above geometries and k (k) values were obtained. Then, the difference of any two k (k) values was calculated. Lastly, all the difference values of k (k) were averaged to obtain the MAE of k Ω 1 − k Ω 2 and k λ 1 − k λ 2 for different ground cover types (shown in Figure 9). λ 1 and λ 2 was taken as the red and near-infrared channel, respectively. A smaller difference means more robustness on geometry. The mean absolute error of k Ω 1 − k Ω 2 is 0.015 while that of k λ 1 − k λ 2 0.1. k Ω 1 − k Ω 2 is much smaller than k λ 1 − k λ 2 . Verification with the RPV model demonstrates the improvement of Equation (8) when applied to the red and near-infrared channels. Furthermore, by addressing the unignorable BRDF issue through the k ratio, the AOD retrievals could be improved not only because the more reasonable assumption of the unchanged k ratio during a certain period of time than the tradition ratio k, but also because the AOD retrieval using the k ratio is more robust to the error than the conventional algorithm ( Figure 6).

Determination of the Surface Reflectance Ratio k
After the AVHRR data preprocess, the clear day was searched by pixel from ±15 days of the retrieval image with the assumptions that the surface bidirectional properties do not change in 31 days and that at least one clear day exists in the 31-day span. The 31-day span was selected because too long a time period may have a significant surface reflectance change. For the contrasting too short a time period there might be no enough cloud clear days to select and the retrieved AOD values may be lower than the ground measurements [22,52]. The retrieval region was confined to non-bright areas where NDVI is greater than 0.15. Then, the clear day was determined by the least TOA reflectance ratio of 0.63 to 0.85 µm because the atmospheric effect would increase the TOA reflectance at 0.63 µm and decrease the TOA reflectance at 0.85 µm in the non-bright areas.
After the clear day was found, the surface reflectance at 0.63 and 0.85 µm of the clear day can be acquired via removing the background AOD and Rayleigh scattering. Knapp et al. [53] demonstrated that the background AOD value error can be introduced to AOD retrieval, but the error is not enlarged in AOD retrievals. In this study, daily Aqua MODIS Deep Blue AOD product (MYD08_D3 Version 6.1) was averaged to produce monthly background AOD values, similar to Hsu et al. [22]. Then, the surface reflectance ratio k of the clear day pixel was determined according to Equation (4). Assuming the surface bidirectional reflectance properties do not change during a period of time, the surface ratio k of the retrieval day can be derived according to Equation (8).

The Dual-Channel AOD Retrieval Algorithm
The flowchart of the AOD retrieval algorithm after the determination of the surface reflectance ratio k is shown in Figure 10. The AVHRR AOD value was retrieved in the 10 × 10 pixels aggregation box (1.1 degree × 1.1 degree). In the box, aerosol model was treated as unchanged [54] assuming that aerosol optical properties change slowly in space [55]. The six aerosol models in Eastern Asia clustered from ground-based measurements [56] were used. Then for every aerosol model in each box, the AOD spatial distribution was retrieved by the following steps. Firstly, a group of 21 AOD initial values (0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, and 5) were constructed. Then for each AOD initial value, the surface reflectance value at 0.63 and 0.85 µm was retrieved by linear interpolation using the 6SV-calculated Lookup Table  (LUT; the 6SV software cf. [41]). In the LUT, the TOA reflectance at 0.63 and 0.85 µm was calculated offline as a function of solar zenith angle, satellite zenith angle, relative azimuth angle, AOD, surface reflectance, and aerosol model. Detailed description of the LUT is shown in Table 3. Therefore, the group of surface reflectance ratio k (having 21 values corresponding to the AOD initial values) was obtained. Figure 10. Flowchart of the dual-channel AOD retrieval algorithm after the determination of the surface reflectance ratio k . It followed that the k group was compared with the surface reflectance ratio k of the retrieval day determined from Section 3.3 to obtain AOD value at 0.55 µm by linear interpolation. The basis for this step is that k would be monotonically decreasing as the AOD value increases. It is because for given surface reflectance values at 0.63 and 0.85 µm, the atmospheric effect would increase the TOA reflectance at 0.63 µm and decrease the TOA reflectance at 0.85 µm in the non-bright areas. Therefore, when the TOA reflectance values are given, the surface reflectance value at 0.63 µm will decrease with the increasing AOD values, while the surface reflectance value at 0.85 µm will increase. Thus, the ratio k between 0.63 and 0.85 µm would be monotonically decreasing as the AOD value increases. If we draw a horizonal line near around 0.8 in the Figure 5, it is also shown that when the TOA reflectance ratio is given (~0.8), the AOD value increases with the decreasing surface reflectance ratio (k ) value (from 0.845 to 0.1177).
Lastly, the final AOD value at 0.55 µm and surface reflectance value at 0.63 and 0.85 µm of every box were retrieved according to the valid result count in the box and the least difference between the retrieved and input surface reflectance ratios. It is noted that the spatial resolution of retrieved AOD and surface reflectance is 11 km.

Uncertainty Study
AOD retrievals with different added errors of k were simulated in this part. Errors ranging from −0.05 to 0.05 stepped by 0.01 were added to the surface ratio k . The value of 0.05 was used because the differences of k for most ground cover types are less than 0.05 as shown in Figure 9. The sun-view angles were the same as Figures 4 and 5. The AOD value was set as 0.4 and 1.5, respectively. A typical vegetation example is shown (the surface reflectance at 0.63and 0.85 µm is 0.05 and 0.46, respectively) in Figure 11. It is revealed that the surface reflectance value at 0.85 µm decreases slowly, whereas the AOD value decreases as the surface reflectance value at 0.63 µm increases. It should be noted that the surface reflectance input at 0.63 µm is much less than that at 0.85 µm. It implies that the k error has a larger effect on the surface reflectance at 0.63 µm than on the surface reflectance at 0.85 µm. When the AOD value is 0.4 and 1.5, the maximum relative errors of AOD are 71.4 and 21.3%, respectively (The relative error δ τ of AOD is defined as here τ ref represents the reference value of τ and as to Figure 11a,b τ ref is 0.4 and 1.5, respectively). Therefore, it is indicated that the accuracy of k has a greater impact on the AOD retrieval when the AOD value is low.

Results and Validation
First, the retrieved surface reflectance value was compared with the atmospherically corrected surface reflectance value using AERONET observation and MOD09 product. Then, the retrieved AOD value was compared with AERONET observation and MODIS AOD value. The spatial distribution of one day is also shown.

Comparison of the Retrieved Surface Reflectance against Atmospherically Corrected Surface Reflectance and MOD09 Product
The retrieved surface reflectance value was compared with the AVHRR atmosphericcorrected surface reflectance value first. The AVHRR atmospheric-corrected surface reflectance was derived with the input of the ground-based AERONET AOD, size distribution, and complex refractive index product and can be considered as benchmark surface reflectance. Figure 12 is the comparison of the surface reflectance at 0.63 and 0.85 µm. The comparison collocations of 0.63 µm shows a more dispersed distribution (R = 0.59) than that of 0.85 µm. It could be caused by the estimated k error, confirming that the k error has a larger effect on the surface reflectance at 0.63 µm than on the surface reflectance at 0.85 µm. Furthermore, the undetected cloud and cloud shadow would intensify this error. The root mean square error (RMSE) of 0.63 and 0.85 µm is 0.03 and 0.029, respectively, also indicating the surface reflectance retrieval at 0.85 µm is more robust than that at 0.63 µm which is consistent with the uncertainty conclusion in Section 3.5. satellite. The channel 2 result was not used for comparison because the spectral response curve of channel 2 of AVHRR is much wider than that of MODIS ( Figure 2). The MOD09GA result which is identified with cloud pollution was removed. The correlation coefficient of channel 1 is 0.76, and the RMSE error 0.031 ( Figure 13). The AVHRR channel 1 surface reflectance is generally less than the MODIS surface reflectance. The difference would be caused by the different spectral response functions between the AVHRR and MODIS bands ( Figure 2). The comparison between AVHRR and MOD09GA surface reflectance results shows that the alignment between the two surface reflectance results is relatively good.  Next, the AVHRR surface reflectance retrieval was compared with daily Terra and Aqua MODIS MOD09GA surface reflectance product of a 500 m spatial resolution in the AERONET stations. The AVHRR channel 1 result of the morning satellite (NM) was compared with the MOD09GA channel 1 result of Terra satellite, while the AVHRR result of the afternoon satellites (NN and NP) was compared with the MYD09GA result of Aqua satellite. The channel 2 result was not used for comparison because the spectral response curve of channel 2 of AVHRR is much wider than that of MODIS (Figure 2). The MOD09GA result which is identified with cloud pollution was removed. The correlation coefficient R of channel 1 is 0.76, and the RMSE error 0.031 ( Figure 13). The AVHRR channel 1 surface reflectance is generally less than the MODIS surface reflectance. The difference would be caused by the different spectral response functions between the AVHRR and MODIS bands ( Figure 2). The comparison between AVHRR and MOD09GA surface reflectance results shows that the alignment between the two surface reflectance results is relatively good.
surface reflectance results shows that the alignment between the two surface re results is relatively good.   Figure 14 shows a comparison map example of AVHRR surface reflectance at 0.63 μm compared with the MYD09 surface reflectance product at 0.66 μ September 2009. Note that the AVHRR algorithm is applicable to non-bright reg the exception of oceans. The scatter-density plot shows similar performance betw AVHRR and MYD09 surface reflectance products with the 0.81 and the RM 0.034 ( Figure 15). The lower AVHRR surface reflectance values than MODIS woul to the different spectral response functions between the AVHRR and MODIS ban  Figure 14 shows a comparison map example of AVHRR surface reflectance retrieval at 0.63 µm compared with the MYD09 surface reflectance product at 0.66 µm on 24 September 2009. Note that the AVHRR algorithm is applicable to non-bright region with the exception of oceans. The scatter-density plot shows similar performance between the AVHRR and MYD09 surface reflectance products with the R 0.81 and the RMSE error 0.034 ( Figure 15). The lower AVHRR surface reflectance values than MODIS would be due to the different spectral response functions between the AVHRR and MODIS bands.

Comparison of the AOD Retrieval against the AERONET and MODIS Prod
We compared the derived AVHRR AOD result with the AERONET data and product taking the spatiotemporal effect into consideration. A mean AE observation within ±30 min of the satellite overpass was calculated so that it c compared with the mean of 5 × 5 pixels centered over the AERONET station [ result is considered as valid if half of the pixels have values. Figure

Comparison of the AOD Retrieval against the AERONET and MODIS Products
We compared the derived AVHRR AOD result with the AERONET data and MODIS product taking the spatiotemporal effect into consideration. A mean AERONET observation within ±30 min of the satellite overpass was calculated so that it could be compared with the mean of 5 × 5 pixels centered over the AERONET station [44]. The result is considered as valid if half of the pixels have values. Figure 16 depicts the relationship between the AVHRR AOD value and AERONET AOD observation at 0.55 µm, yielding the linear regression equation of y = 0.074 + 0.89x with an R of 0.88, an RMSE of 0.15, and 83% of collocations within error lines (±0.15 ± 0.15τ, τ is AOD). The non-zero intercept may be attributed to the calibration inaccuracies or the surface tackling error, and the deviation of the slopes from unity results from the inappropriate selection of aerosol model in the AVHRR algorithm [44]. The RMSE error is probably due to the random error caused by the surface reflectance changes or subpixel cloud contamination, while the R denotes the possible ability to detect aerosol signal from the AVHRR sensor [22,[57][58][59][60].
Remote Sens. 2021, 13, x FOR PEER REVIEW Time series comparison of AVHRR AOD retrievals versus AERONET measurements within ±30 min of the satellite overpass is shown in Figure 17. Th AVHRR onboard satellites were involved, namely NM, NN, and NP. It is shown t AVHRR AOD retrievals can capture the AOD changes over time. Time series comparison of AVHRR AOD retrievals versus AERONET AOD measurements within ±30 min of the satellite overpass is shown in Figure 17. The three AVHRR onboard satellites were involved, namely NM, NN, and NP. It is shown that the AVHRR AOD retrievals can capture the AOD changes over time. Time series comparison of AVHRR AOD retrievals versus AERONET AOD measurements within ±30 min of the satellite overpass is shown in Figure 17. The three AVHRR onboard satellites were involved, namely NM, NN, and NP. It is shown that the AVHRR AOD retrievals can capture the AOD changes over time. The conventional estimate of AOD at 0.55 µm using 0.63 µm channel was also retrieved. The conventional estimate of AOD used the same clear day results as the k ratio. Validation against the AEROENT measurement in Figure 18 shows that the linear regression equation is y = 0.17 + 0.57x with an R of 0.78 and an RMSE error of 0.22, which is inferior to the AOD retrieval using the k ratio. One example of the large-scale AOD retrievals from AVHRR and MODIS is shown in Figure 19.  The text describes the linear regression relationship, R, RMSE error, and ratio within the error lines.
One example of the large-scale AOD retrievals from AVHRR and MODIS is shown in Figure 19. Both the AQUA/MODIS Dark Target and Deep Blue AOD products are utilized. The AVHRR AOD spatial distribution is similar to the MYD04 product, which is low in 42-54 •

Conclusions
The orbit drift and phasing running of NOAA satellites and large view geometry variation would cause an unignorable BRDF effect in the determination of surface reflectance from time series of AVHRR measurements. Our forward sensitivity study found the need for dual channels at 0.63 µm and 0.85 µm to retrieve AOD instead of the singlechannel algorithm using 0.63 µm. We used the surface reflectance ratio k between 0.63 µm and 0.85 µm to address the surface BRDF effect. The ratio was demonstrated from both the forward-sensitivity and theoretical point of view and the validity was examined. It was found that the AOD inversion using the k ratio is more robust to the input error than the conventional 0.63 µm-channel algorithm. Moreover, it is more reasonable to assume that the k ratio remains unchanged during a certain period of time compared to the traditional ratio k when addressing the BRDF issue. Then, the AVHRR dual-channel algorithm was proposed to determine AOD at 0.55 µm and surface reflectance at 0.63 and 0.85 µm. This novel approach was applied to Eastern Asia. The surface reflectance retrieval was compared with the atmospherically corrected surface reflectance value and MODIS surface reflectance product. The surface reflectance at 0.85 µm has better statistic results than that at 0.63 µm, consistent with the uncertainty conclusion. The AOD retrieval was compared with AERONET AOD observation, yielding the linear regression equation of y = 0.074 + 0.89x with an R of 0.88, an RMSE error of 0.15. The fraction of the AOD validation collocations within the error lines (±0.15 ± 0.15τ) was 83%. The favorable correlation provides confidence in the utility of our AVHRR AOD algorithm over land. Time series comparison between AVHRR AOD retrievals and AERONET AOD measurements shows that the AVHRR AOD retrievals can capture the AOD changes over time. The proposed dual-channel AOD algorithm taking into account the surface BRDF effect is proved outperforming the conventional 0.63 µm-channel method.
The generic nature of the dual-channel algorithm makes it appropriate for AOD retrievals from common sensors, especially from those that lack the shortwave infrared channel that is required in the MODIS Dark Target AOD algorithm. Nonetheless, the present algorithm is confined to non-bright regions because of the use of the clear day. Moreover, the algorithm is more suitable for the situation in which the surface reflectance does not change dramatically. Further research could focus on improving the gathering of the surface reflectance of an adjacent clear day, such as adopting the quantitative expression of surface bidirectional reflectance [46] and employing a non-Lambertian forward model [61,62].
AVHRR may provide an almost 40-year unique archive of AOD retrieval over land from the 1980s to present, during which time aerosol plays a particularly important and unclear role in the global and regional climate. To achieve this objective, more studies are still needed; for example, earlier aerosol model assumption and AOD validation without the AERONET observation available would be helpful.