The Classification of Reflectance Anisotropy and Its Application in Albedo Retrieval

: The land surface albedo reflects the ability of the surface to reflect solar radiation and is a critical physical variable in the study of the Earth’s energy budget and global climate change. Algo ‐ rithms for the retrieval of albedo usually require multi ‐ angle measurements due to surface anisot ‐ ropy. However, most of the satellites cannot currently provide sufficient and well ‐ distributed ob ‐ servations; therefore, the accuracy of remotely sensed albedo is constrained. Based on the Moderate ‐ Resolution Imaging Spectroradiometer (MODIS) Bidirectional Reflectance Distribution Function (BRDF) and albedo product (MCD43A1), this study proposed a method to further subdivide reflec ‐ tance anisotropy and build an updated database of BRDF archetype, using both the Anisotropic Flat Index (AFX) and Perpendicular Anisotropic Flat Index (PAFX). The BRDF archetypes were used to fit the corresponding MODIS BRDF, and the optimal number of BRDF archetype categories was determined according to the tendency of fitting error. The effect of surface anisotropy and observa ‐ tion noise on albedo retrieval were explored based on simulated MODIS reflectance. Finally, the BRDF archetype A2P2 was taken as prior knowledge to retrieve albedo from a different number of MODIS observations, and the result was validated by the high ‐ quality MODIS albedo product. The results show that the fitting error between BRDF archetypes and MODIS BRDF shows a rapid de ‐ cline when introducing the PAFX in the classification process. A 3 ‐ by ‐ 3 matrix of BRDF archetypes, which occupy 73.44% and 70.13% of the total decline in the red and NIR band, can be used to rep ‐ resent the characteristics of reflectance anisotropy. The archetype A2P2 may be used as prior knowledge to improve the albedo retrieval from insufficient observations. The validation results based on MODIS observations show that the archetype A2P2 ‐ based albedo can reach root ‐ mean ‐ square errors (RMSEs) of no more than 0.03.


Introduction
Albedo is defined as the ratio of the reflected solar ration to the total incoming solar radiation over the whole solar spectrum and is a key factor in the study of global environmental change [1,2]. The reflectance characteristics of most of the surface types are anisotropic. The surface reflectance anisotropy, which is usually described by the Bidirectional Reflectance Distribution Function (BRDF) [3,4], is dependent on the surface structure and the angles of incidence and reflection. Albedo is the integral of reflectance over the incident and emitted hemispheric space, so accurate albedo retrieval needs to consider the anisotropic characteristics of the surface [5,6].
A sufficient number of multi-angular observations describing the anisotropy characteristics of the surface is usually required to estimate the albedo with good quality [7].
Currently, the most widely used satellite multi-angular observations are from the MOIDS and Polarization and Directionality of the Earth's Reflectance (POLDER), which have large-view zenith angles along and across the track. The multi-angular observations from these sensors can only be obtained by the accumulation of sequential observations over a specified period [7,8]. Based on the kernel-driven BRDF model, the MODIS BRDF and albedo product are retrieved [3,9] and widely used in studying surface reflectance anisotropy and surface albedo change [10][11][12][13]. However, for medium-and high-resolution satellites, it would be difficult to obtain the multi-angle observations due to the small-view zenith angles, the long revisit period, and the influence of clouds.
Prior BRDF knowledge may be able to improve the accuracy of albedo retrieval from insufficient observations, or even one single directional observation [14]. Varying prior BRDF knowledge has been used in albedo retrieval [10,15,16] and reflectance normalization [17]. Based on the prior BRDF of different land cover classes, Strugnell and Lucht [18] presented a backup BRDF/albedo algorithm when the number of observations is not sufficient. Based on the Bayesian theory of inference, a priori knowledge was introduced to effectively solve the problem of retrieval instability [14]. Associating the reflectance anisotropy to land cover types or normalized difference vegetation index (NDVI), different methods have been proposed to retrieve 30 m Landsat surface albedo according to the prior BRDF extracted from MODIS [10,15,19]. Landsat albedo was also retrieved based on direct estimation algorithm [20,21]. For the study of reflectance normalization, reflectance data have been normalized to the nadir based on the surface reflection anisotropy information extracted through POLDER or MODIS data [17,22,23]. To classify the complex reflectance anisotropy, the BRDF model parameters were normalized, and a BRDF archetype database was established based on the Anisotropic Flat Index (AFX) and MODIS BRDF product [24]. The accuracy of the BRDF archetypes to represent the reflectance anisotropy in albedo retrieval from multi-angular reflectance has been discussed [25].
Extracting effective information from large amounts of historical BRDF products can improve the understanding of surface reflectance anisotropy and its effect on albedo retrieval. Based on the AFX theory, a new Perpendicular Anisotropic Flat Index (PAFX) was proposed in this study. These two indices were combined to classify the surface anisotropy characteristics, and a new BRDF archetype matrix was built. The BRDF archetypes were used as prior knowledge to explore the effect of reflectance anisotropy on albedo and the possibility of using one specific BRDF archetype for albedo retrieval.

Kernel-Driven BRDF Model
The kernel-driven BRDF model was generally used to retrieve the surface albedo from multi-angle measurements, and a linear combination of kernels of a certain physical meaning was used to describe the bi-directional reflection properties of the surface [26][27][28]. The model takes the form of Equation (1): where fiso(λ), fvol(λ), and fgeo(λ) are the spectrally dependent model parameters, indicating the percentage of isotropic scattering, geometric-optical scattering, and volumetric scattering components in the reflectance, respectively. Kvol is derived from volume scattering radiative models [29], and Kgeo is retrieved from surface scattering and geometric shadow casting theory [30], and they are functions of the incident and observed angles. R (θi, θr, ψ, λ) is a bidirectional reflectance distribution function in waveband λ, which is the illumination zenith angle, θi, the view zenith angle, θr, and the reflectance at the relative azimuth angle, ψ.
The BRDF model is a combination of the volumetric scattering kernel and the geometric-optical scattering kernel, and the different kernel combinations produce different results. The kernels chosen in this study are the Ross-Thick and the LiSparseR. It has been proved by many studies that such a combination has the best fitting effect, as well as extrapolation ability to estimate the surface reflectance with high accuracy and convenience [3,26,28]. The Ross-Thick kernel has an improved effect on the reflectivity when the view zenith angle is large and the BRDF shape is "bowl-shaped", while the LiSparseR kernels represent a somewhat dome-shaped BRDF curve in the backward-scattering direction [26]. The actual BRDF shape is expressed by these two shapes after adding weights, so the different BRDF shapes can indicate the characteristics of reflectance anisotropy. The BRDF model parameter (fiso, fvol, and fgeo) can be generated by fitting the multi-angular observation with the Ross-Thick/LiSparseR model by linear regression. Then the bi-directional reflection in any viewing and illumination directions can be obtained through the model parameters and the kernels in the forward mode.

MODIS Database
MODIS surface reflectance products MOD09GA and MYD09GA are from Terra and Aqua satellites, respectively. The product provides surface reflectance data from nearsurface measurements, independent of atmospheric absorption and scattering [6]. Based on a combination of the MODIS operational algorithm and 16 days of atmospherically corrected surface reflectance data, the MODIS product provides a global 500 m resolution BRDF/albedo product with a long time series [3]. In particular, MCD43A1 includes BRDF model parameter products for each band (with parameters fiso, fvol, and fgeo), MCD43A2 describes its quality, and MCD43A3 includes the black sky albedo (BSA) and white-sky albedo (WSA) products. Many studies have shown that the MODIS BRDF/albedo products can accurately reflect the conditions of the surface and meet the precision requirements in different situations [13,31,32].
The global time series MODIS BRDF parameter and the quality products in 2015 in the red and near-infrared (NIR) bands were used in this study. The global equal interval sampling method was applied to the MODIS data, with an interval of 2°. The number of selected pixels is about 15,000, and the number of high-quality time-series MODIS BRDF parameter samples is about 2 million. The samples contain most of the IGPB land-cover types and include a wide variety of surface anisotropy reflectance characteristics.
To further verify the applicability of the extracted BRDF archetype in different observations, MODIS surface reflectance (MOD09GA) at two periods (2021.101-116 and 305-320) in tile h20v11 was also used in this study. This tile is located in the southern part of Africa, with the main surface types being grassland, open shrubland, savanna, and cropland/natural vegetation mosaic. Only the first layer of the surface reflectance with high quality during the 16-day period was used. Figure 1 shows an example of the MODIS angular samplings at the two periods. The observations during the first period are close to the cross-principal plane (CPP), with the solar zenith angle being around 45°, while the observations during the second period are located near the principal plane (PP), with a solar zenith angle range of 13° to 32°. The surface albedo was retrieved from the MODIS surface reflectance and the BRDF archetypes, and the accuracy of albedo was validated by the high-quality MODIS albedo product (MCD43A3).

AFX and PAFX
The extraction of prior BRDF knowledge requires the classification, effective extraction, and evaluation of information about surface reflectance anisotropy characteristics, with the help of some classification indices. AFX is an angular index constructed for the anisotropic reflection pattern of the MODIS sensor and is defined as the ratio of the WSA to the surface isotropic parameter, fiso [24].
The WSA, which can be obtained by integrating the reflectance over both the viewing hemisphere and illumination hemisphere, can be derived as Equation (2), according to the kernel-driven BRDF model. The AFX that relates to the geometric-optical scattering and volumetric scattering weights can be written as Equation (3).

,
( where Kv and Kg are the directional-hemispherical integrals of the BRDF model kernels, and the values are 0.189184 and −1.137762, respectively. Combining Equation (3) and the shapes of the Ross-Thick kernel and LiSparseR kernel in the PP, it can be deduced that different values of AFX can show various BRDF shapes. When the volumetric scattering effect is larger than the geometric-scattering effect (AFX > 1), the shape of BRDF shows a bowl shape with a low center and high sides. When the volumetric scattering effect gradually decreases and the geometric-optical effect occupies the dominant position (AFX < 1), the shape of the BRDF approximates a dome shape. When the two effects are close to each other, the AFX is close to 1, and the shape of BRDF is flat [24]. AFX, which is well-associated with the shape of BRDF, provides a new method for recognizing complex anisotropy of the surface and can be used to quantitatively evaluate the characteristics of reflectance anisotropy.
In the previous study, the BRDF archetypes were extracted from normalized BRDF model parameters, which removed the effect of spectral reflectance amplitude and can be compared directly [24]. The BRDF model parameters were normalized by a scale factor of 2×fiso to limit the range of normalized BRDF shapes to 0-1. The normalized BRDF model parameters (F) can be written as follows: Then Equation (3) can be rewritten as follows: Take the Fvol and Fgeo as independent and dependent variables, and an F parameter space could be constructed. Setting a specific value of AFX, the parameter space would be divided into two parts (white and gray parts in Figure 2). In Figure 2, the black line has the same AFX value; the AFX value in the upper part (white part) is smaller; and in the bottom part (gray part), the AFX value is bigger. The variation in the BRDF shapes in one archetype class defined by AFX is still notable due to the large range of model parameters.
To subclassify the surface reflectance anisotropy, a new index PAFX was introduced to redefine the archetype class. PAFX, which is perpendicular to Equation (5), can be written as Equation (6). Similarly, given a specific value of PAFX, the parameter space would be divided into two parts (left part and right part in Figure 2). When combining the two classification results, the parameter space would be divided into four classes (A1P1, A1P2, A2P1, and A2P2) in this condition. The variation in the model parameter in the new category is much smaller than that in the original category. The optimal class number for AFX and PAFX categories are discussed later in the paper. 2 2 (6) Figure 2. The relationship between the two indices in parameter space. The region name rule is the category index number when each of AFX and PAFX is classified independently; for example, A1P1 means that the region belongs to the first category of AFX and the first category of PAFX.

Classification Strategy and Evaluate
Previous research has shown that the AFX can indicate the pattern of reflectance anisotropy, which can be used to distinguish surface reflectance anisotropy and benefit the prior BRDF extraction [24]. The extraction of surface anisotropy information by AFX mainly uses the Iterative Self-Organizing Data Analysis Techniques Algorithm (ISO-DATA) that generates n categories of classification schemes based on the numerical classification of AFX. The mean values of the model parameters within a category were used as the model parameters of the BRDF archetype. The number of classes is determined by the difference between archetype BRDF and MODIS BRDF. The classification method used in this study is similar, and the main difference is that both AFX and PAFX are used in the classification process. AFX and PAFX are independently used as classification indicators for BRDF characteristics. The combination of the two classification results establishes the final classification result. The parameter space would be divided into an nby-m matrix of categories, and the mean value of the model parameter in each category is taken as parameter for the BRDF archetypes.
The BRDF archetype (BRDF') is used to fit the corresponding MODIS BRDF, according to the least-squares method, to evaluate the classification. In this process, the simulated MODIS multi-angular reflectance (ρ), which has 1440 observations (the solar zenith angle ranges from 0° to 70°, with an interval of 5°; the view zenith angles range from 0° to 70°, with an interval of 10°; and the relative azimuth angles range from 0° to 360°, with an interval of 30°), is used. First, the reflectance (ρ') that has the same observation geometry with ρ was simulated from the kernel-driven BRDF model and the BRDF archetypes: Then we used an adjustment, a, to minimize the difference between MODIS BRDF and adjusted BRDF archetypes, using the technique of least-mean squares: The adjusted BRDF archetypes and the fitting error (RMSEa) of each sample can be written as follows: Finally, the mean value of RMSEa is used to evaluate the difference between BRDF archetypes and MODIS BRDF and the representation of BRDF archetypes to various reflectance anisotropy. Beyond that, the WSA can be calculated through the adjusted BRDF archetypes and Equation (3).

WSA Retrieval from BRDF Archetype
In this study, the adjusted BRDF archetypes that can best fit the observations can be calculated through Equations (7)- (9). The WSA of the BRDF archetype (WSA') can be calculated by using Equation (2), and the BRDF archetype-based WSA can be written as follows: , The high-quality MODIS BRDF product was used to evaluate the accuracy of BRDF archetype-based WSA (WSAb). The RMSEr and BIAS between the high-quality MODIS WSA (WSAa) and BRDF archetype-based WSA (WSAb) is used to describe the accuracy of the BRDF archetype. , where n is the number of the samples.

Effect of Reflectance Anisotropy and Observation Noise
The BRDF archetypes cover all the reflectance anisotropy characteristics, which can facilitate the study on the impact of reflectance anisotropy on albedo retrieval from remotely sensed observations. The number of observations and their distribution characters are usually the main factors considered in albedo retrieval from observations [7]. Based on this analysis, we designed three different types of observations based on simulated MODIS reflectance. The first type is sufficient multi-angular reflectance, which is distributed in PP and CPP, with an interval of 10° and the maximum view zenith angle of 60° ( Figure 3). The second type is the insufficient multi-angular reflectance, which has three observations in the PP or CPP, with the largest view zenith angle of 15°. The third type is the only one single directional reflectance in any direction of the viewing hemisphere, with the largest view zenith angle of 60°. The solar zenith angle is set to 30° for the three types. Different BRDF archetypes are used to fit the simulated reflectance to generate albedo (Equations (7)-(11)). The retrieved albedo is compared with the high-quality MODIS albedo product to show the effect of reflectance anisotropy and to evaluate the accuracy of albedo. Since the real observations contain a certain degree of noise, which is not considered in the simulated reflectance, it is necessary to add random noise (the largest proportion is 20%) to the simulated reflectance to learn the effect of noise on albedo retrieval and further illustrate the feasibility of using BRDF archetypes for albedo retrieval. In all of the above cases, the RMSEr and BIAS between MODIS albedo and BRDF archetype-based albedo are used to quantitatively evaluate the accuracy of albedo.

BRDF Archetype Database Based on AFX/PAFX
The AFX and PAFX values for each sample can be calculated from high-quality MODIS BRDF parameter products (MCD43A1). To determine the best number of classifications, first, the AFX and PAFX values are classified into 1 to 10 classes, using ISODATA, according to the previous study [24]. Then the two classification results are combined to reclassify the data. The BRDF archetype of each new class can be easily calculated. Finally, the BRDF archetypes are used to fit the MODIS BRDF within the class. The mean fitting RMSE (RMSEa) of all samples is used to measure the representativeness of BRDF archetypes to reflectance anisotropy.
The combination of AFX and PAFX can generate a matrix RMSEa sized 10 by 10. Figure 4 shows the change in RMSEa with the number of classifications. It shows that, as the number of AFX and PAFX classifications increases, the RMSEa decreases gradually. As compared with the RMSE only based on AFX, the decrease in RMSEa is more pronounced. The RMSEa showed a marked decline with the increase in classification number at the beginning, while the declining tendency becomes slow when the number of classifications becomes large. The results suggest that the combination of AFX and PAFX can effectively improve the classification accuracy of reflectance anisotropy.  Considering the quality of classification results and the convenience of BRDF archetype application, we chose the combination of three classification results by AFX and PAFX as the final classification results. At this time, the RMSEa is 0.0092 and 0.0142, which occupy 73.44% and 70.13% of the total decline (the maximum minus the minimum) in the red and NIR band. The 3-by-3 matrix of BRDF archetypes can better present the surface reflectance anisotropy characteristics. The ranges of AFX and PAFX for different classes in the red and NIR bands are shown in Table 1.   Figure 6 shows the three-dimensional BRDF shapes of the 3-by-3 matrix of BRDF archetypes in red and NIR bands. Each band contains nine BRDF archetypes representing different surface anisotropy characteristics. Table 2 shows the model parameters for each BRDF archetype in Figure 6. The A-series BRDF archetypes are distinguished by AFX proposed by previous research [24]; the pattern of BRDF shape gradually changes from roofshaped, with mainly geometric-optical scattering, to bowl-shaped, with mainly volumetric scattering, with the class number change from A1 to A3.
AFX can effectively indicate the main anisotropic reflectance patterns at the Earth's surface, while the PAFX proposed in this study has the ability to detail the BRDF shapes within an A-series class. With the class number change from P1 to P3, both of the two parameters increase gradually. The BRDF shapes of P3 have a relatively large (small) value in the backward (forward) directions as compared with P1. The archetype A1P3 represents the strongest reflectance anisotropy, with the BRDF shapes having the largest difference, especially in the forward and backward directions. The archetype A3P1 is flatter since the model parameters are close to 0, and this archetype can be considered to be the Lambertian assumption. The archetype A3P3 has an obvious bowl shape, with the strongest volumetric scattering and a relatively small geometric scattering. The BRDF archetypes in the red and NIR bands have high consistency, and this is because the effect of spectral information is removed in the BRDF archetypes, which only keep the pure reflectance anisotropy characteristics.

Albedo Retrieval from a Specific BRDF Archetype and Multi-Angular Reflectance
Based on Section 2.4.2, we learned the effect of surface reflectance anisotropy and observation noise on albedo retrieval for different situations by using the extracted BRDF archetypes. Different BRDF archetypes were used to fit the simulated multi-angular reflectance, using the least-squares method. Figure 7 shows the effect of different surface anisotropies on albedo retrieval from sufficient and insufficient multi-angular reflectance with a solar zenith of 30°.
As compared with the MODIS albedo product, the WSA retrieved from A1-series archetypes usually underestimate the albedo, especially for the archetype A1P1, while the WSA retrieved from A3-series archetypes usually overestimates the albedo, especially for the archetype A3P3. The WSA retrieved from A2-series archetypes has a high consistency with the MODIS albedo product. Referring to the shapes of BRDF archetypes in Figure 6, this pheromone could mainly be caused by the difference in BRDF shapes in the forward directions. The albedos retrieved from archetype A2P2 are more stable than the other archetypes, with the smallest RMSEr and a BIAS nearest to 0 in all three multi-angle observation situations.
For the insufficient multi-angular reflectance in CPP, similar results also can also be obtained with a relatively larger RMSEr and BIAS as compared with sufficient reflectance. The retrieved albedo from the insufficient reflectance in CPP has a similar performance with the reflectance in PP, and the results are not shown.

Albedo Retrieval from a Specific BRDF Archetype and a Single Directional Reflectance
The effect of reflectance anisotropy on albedo retrieval from a single directional reflectance is discussed in this section. Figure 8 shows the comparison of MODIS albedo with albedo retrieval from nine BRDF archetypes and nadir reflectance in the red and NIR bands. Similar to the results of multi-angular observations, the A1-series archetypes would underestimate the albedo, while the A3-series archetypes would overestimate the albedo. It also shows that the WSA retrieved from archetype A2P2 has a good agreement with MODIS albedo, with an RMSEr of 0.021 and 0.036 and a BIAS of 0.002 and 0.006 in the red and NIR bands, respectively. The retrieved albedo from A2P2 is slightly less accurate than the albedo retrieved from multi-angular reflectance, but the result still has good consistency with MODIS albedo.
For different observation angles, the effect of the reflectance anisotropy on albedo retrieval may be different; we further extended the directional reflectance from nadir to the entire viewing hemisphere. Similarly, the RMSEr in each direction is used to evaluate the possibility of using the BRDF archetypes and the reflectance in the direction. Figure 9 shows the distribution of RMSEr in the viewing hemisphere with a solar zenith angle of 30°, and the small value of RMSEr corresponds to a higher accuracy of albedo. Overall, the archetypes A2P1, A2P2, A3P1, and A3P2 have higher retrieval accuracy, with the RMSEr having a small value (less than 0.05) in the entire viewing hemisphere, while other BRDF archetypes have better retrieval results only in some limited directions. Noteworthily, the albedo retrieved from reflectance on the CPP is relatively stable, while the albedo retrieved from reflectance on the PP is less accurate, especially in the backward and around the hotspot directions. This may be caused by the remarkable changes in reflectance anisotropy on the PP. In the actual situation, observations obtained by satellites may contain a certain level of noise caused by internal and external factors (e.g., instrument senescence, climate, or underlying surface changes). It is necessary to consider the impact of noise when learning the effect of reflectance anisotropy on albedo retrieval. The impact of random noise in the multi-angular observation on albedo could be offset by using the least-squares method. However, when retrieving albedo from BRDF archetypes and a single directional reflectance, the noise in the reflectance would be delivered to the albedo completely, using the proposed algorithm. Figure 10 shows the results of the albedo retrieval from different BRDF archetypes and nadir reflectance with random noise (the largest proportion is 20%).
For the albedo retrieval form nadir reflectance with noise, the A1-series archetypes would underestimate the albedo, while A3-series archetypes would overestimate albedo. The scatter plot is more discrete, and the RMSEr tends to increase as compared to the previous results. The RMSEr of archetype A2P2 increased from 0.021 to 0.038 and from 0.036 to 0.058 in the red and NIR bands. For the entire viewing hemisphere (Figure 11), although the RMSEr values of different BRDF archetypes tend to increase after adding noise, the albedo retrieval from archetype A2P2 retains an acceptable accuracy. The random noise in the observation could reduce the accuracy of albedo retrieval from BRDF archetypes and directional reflectance, and the archetype A2P2 can be used to remove the effect of reflectance anisotropy.
(a) (b) Figure 10. Comparison of MODIS albedo with albedo retrieval from different BRDF archetypes and nadir reflectance with a random noise range from −20% to 20% in the red (a) and NIR (b) bands.
(a) (b) Figure 11. Same as Figure 9, but for the albedo retrieval from directional reflectance with a random noise range from −20% to 20% in the red (a) and NIR (b) bands.

The Application of the Archetype in MODIS Reflectance
To further illustrate the effect of different surface anisotropy reflectance on albedo retrieval, this section uses real MODIS observations for further validation. For the selected MOD09GA data, only the first layer of reflectance was used in this study. As the observations were obtained over the course of 16 days, the number of valid observations for a pixel range from 0 to 16 after filtering by quality flags. Figure 12 shows the number of observations for each pixel at the two periods in the red and NIR bands. It has been proved that the archetype A2P2 has the best ability to correct the effect of reflectance anisotropy in albedo retrieval based on the simulated reflectance. In this section, the archetype A2P2 is further applied to the MODIS observations, and the retrieved albedo is compared with the high-quality WSA data from the MCD43A3 product. Figures  13a,d and 14a,d show the MODIS albedo product and the albedo retrieval from archetype A2P2 and multi-angular MODIS observations in red and NIR bands within 2021.101-116 and within 305-320. To further verify the possibility of using archetype A2P2 to correct the effect of reflectance anisotropy on albedo from a single directional reflectance, albedo (b and e in Figures 13 and 14) is also retrieved from a single directional reflectance that is selected randomly from the multi-angular observations in this process. The results show that, although the albedo based on a single directional reflectance shows a certain degree of discontinuity, both of the two albedos have high spatial consistency with MODIS albedo products.  Figure 15 shows the comparison of retrieved albedo based on the archetype A2P2 with the high-quality MODOS albedo product. When based on multi-angular observations, the RMSEr is 0.007 and 0.12 for red and NIR bands for the data within 101-116, while the RMSEr is 0.016 and 0.024 for the data within 305-325. The RMSEr is slightly greater for the data within 305-325, and this is mainly caused by the relatively significant reflectance anisotropy characteristics around the PP. In comparison, the RMSEr shows a slight increase when there is only a single directional reflectance. In general, the RMSEr of both multi-angular and a single directional observation at the two periods satisfies the accuracy requirement of albedo applications, especially for the observations close to the CPP. The effect of the number of observations is further explored based on the MODIS observations. Albedo is calculated from archetype A2P2 and the same number of observations for different pixels. To make full use of the data, the pixel will participate in the calculation when the specific number of observations is less than or equal to the number of observations of the pixel. The subset of observations is selected randomly from the multi-angular data of the pixels. This strategy could provide sufficient sample for this study and make the retrieval results more reliable. Figure 16 shows the RMSEr between MODIS albedo and albedo retrieved from different numbers of observations in red and NIR bands at two periods. When there is only a single directional reflectance, the RMSEr reaches its maximum. With the increase in the number of observations, the RMSEr values decrease gradually. The results indicate that increasing the number of observations could improve the accuracy of albedo based on archetype A2P2. The RMSEr is always less than 0.02 and 0.03 in the red and NIR bands but it still meets the accuracy requirement. Overall, archetype A2P2 can be taken as prior BRDF knowledge of reflectance anisotropy in the albedo retrieval from insufficient observations, even from a single directional observation.

Discussion
At present, more and more satellites with high resolution have been launched. Most of the satellites cannot provide sufficient multi-angular observations within a short period due to the narrow viewing field and the long revisiting period. It has been a challenge to explore how to make the albedo retrieval from the insufficient number of observations more accurate [3,10,18]. This study focused on the classification of reflectance anisotropy and explored the effect of reflectance anisotropy on albedo retrieval, and it also explored the possibility of using a specific BRDF archetype to retrieve from insufficient observations.
MODIS BRDF products have been widely used for the extraction of surface reflectance anisotropy, but conventional methods need to introduce other indices [18,33] due to the complexity of reflectance anisotropy. Based on the kernel-driven BRDF model, Jiao et al. proposed a method to link the BRDF shape and AFX, which can indicate the change pattern of surface anisotropy reflectance characteristics [24]. AFX was used to classify MODIS BRDF products, and a BRDF archetype database was built. The BRDF shapes within each BRDF archetype class still show a certain degree of discretion due to the large range of model parameters in each class. To solve this problem, we propose a new indicator which is perpendicular to AFX, named PAFX. The combination of AFX and PAFX can classify the reflectance anisotropy into a more delicate classification. With the increase in the number of classifications, the RMSEa declined quickly to the minimal level acquired merely by AFX. Considering the quality of classification results and the convenience of BRDF archetype application, a 3-by-3 matrix of BRDF archetype database is taken as the final classification results. The nine BRDF archetypes can effectively describe all kinds of reflectance anisotropy.
Based on the simulated reflectance, this study has proved that the geometric-dominated A1-series archetypes usually underestimate the albedo, while the volumetric-dominated A3-series archetypes usually overestimate the albedo. The archetype A2P2 can achieve a relatively accurate albedo from insufficient observations or even from a single directional reflectance. This result is consistent with the previous studies [24,34], but it has a higher accuracy. It is important to point out that the scatter plot of the model parameter usually has two gathering centers ( Figure 5). One of the gathering centers that has a small value of the model parameter is away from the archetype A2P2. This phenomenon may conflict with the method that one specific BRDF archetype can be used for all the pixels to correct the effect of reflectance anisotropy in albedo retrieval. The pixels around the gathering center usually have a small NDVI value and a relatively high reflectance, and the land-cover type for the pixels is usually non-vegetated. If the two gathering centers can be divided, the accuracy of albedo retrieval from a specific BRDF archetype may be further improved.
The archetype A2P2 has the best ability to correct the effect of reflectance anisotropy in albedo retrieval for the pixels with various reflectance anisotropy characteristics. That may be because the albedo is the integration of BRDF over both the viewing and illuminating hemisphere. The integration of BRDF may eliminate the difference of reflectance anisotropy in albedo retrieval. However, if a specific directional reflectance would be used, for example, normalizing the directional reflectance to nadir, the applicability of the one BRDF archetype method needs to be discussed in detail. The varying reflectance anisotropy around the CPP is relatively weak, while around the PP, it is relatively significant, and the observation of most satellites is located around the CPP most of the time. In this case, the BRDF archetype may also be applicable in normalizing the directional reflectance. However, if the observations are around the PP, the BRDF archetype may be less effective.
Although the new BRDF archetype database improved the classification accuracy of reflectance anisotropy, the problem of how to link the BRDF archetypes to the surface remains unresolved. Using the archetype A2P2 for albedo retrieval from insufficient observations is only a special case. The application of BRDF archetypes requires more effect.

Conclusions
In this study, PAFX, which was defined according to the AFX, was added to the classification of reflectance anisotropy. Both AFX and PAFX were combined to classify the reflectance anisotropy represented by global MODIS BRDF products. The mean value of the MODIS parameter within a class was taken as the BRDF archetype, and the BRDF archetype was used to fit the MODIS BRDF within the class. The best number of classifications was determined by the changing tendency of fitting RMSEa. The effect of reflectance anisotropy on albedo retrieval from multi-angular or a single directional reflectance was analyzed based on the simulated reflectance. The archetype A2P2 was taken as prior knowledge for the retrieval of albedo from MODIS observations, and the retrieved albedo was compared with the high-quality MODIS albedo product (MCD43A3). The results show the following: (1) The combination of AFX and PAFX can classify the reflectance anisotropy into a more delicate classification. The difference between BRDF archetypes and MODIS BRDF declines rapidly with the increase in the number of categories. A 3-by-3 matrix of BRDF archetypes, which occupy 73.44% and 70.13% of the total decline in the red and NIR bands, respectively, can describe the main characteristics of surface reflectance anisotropy. (2) The BRDF archetype that has relatively strong geometric-scattering or volumetric-scattering characteristics could underestimate or overestimate albedo, while the archetype A2P2 can be used to improve the accuracy of albedo retrieval from insufficient multi-angular reflectance or even a single directional reflectance. The validation results based on MODIS observations show that the accuracy of albedo retrieval from archetype A2P2 improved continuously with the different number of observations and can achieve RMSEs of no more than 0.03.
The algorithm proposed in this study provides a more delicate reflectance anisotropy classification, and whether the currently extracted BRDF archetype database can be further optimized and the application field of the BRDF archetypes should both be studied in the future research.  Data Availability Statement: All satellite remote sensing and field measured data used in this study are openly and freely available. The Collection 6 MODIS BRDF parameter and land-cover products are available at https://search.earthdata.nasa.gov/search, accessed on 1 September 2021.