## 2 Methodology
### BRDF Archetype Inversion Algorithm  
The core idea of surface albedo inversion using the BRDF archetype is to construct a BRDF archetype database based on a large amount of multiangle reflectance observation data and establish a linear relationship between the BRDF archetype and real BRDF to reduce the amount of input data and provide stable inversion results. The specific steps include calculating the fitting error between the reflectivity corresponding to each BRDF archetype and the real surface reflectivity under the same geometric conditions and selecting the BRDF archetype corresponding to the minimum error as the optimal BRDF archetype. Here, it is assumed that there is a set of multiangle reflectivity observation data \(\rho:\{\rho_{1},\rho_{2},\cdots,\rho_{n}\}\), which corresponds to the \(BRDF\) of the real surface. It is assumed that the corresponding BRDF archetype is known as \(BRDF^{\prime}\), and the reflectivity data corresponding to the observation position of the data \(\rho\) are \(\rho^{\prime}:\{\rho_{1}^{\prime},\rho_{2}^{\prime},\cdots,\rho_{n}^{\prime}\}\). Taking the adjustment coefficients A and B, the difference between \(\rho\) and \(A\cdot\rho^{\prime}+B\) is minimized by translating the BRDF archetype. By calculating the minimum regression fitting error between the real surface reflectance and the BRDF archetype fitting reflectance, the calibration adjustment coefficients A and B are obtained. Since the conventional least squares fitting method is greatly affected by abnormal observations, to ensure the inversion accuracy and improve the inversion stability, the Huber loss fitting method is introduced. This method combines the advantages of root mean square error (MSE) and mean absolute error (MAE) fitting and reduces the sensitivity to outliers. By minimizing the loss function, the optimal fitting coefficient can be obtained. The core formula is  
\[\min_{w,\sigma}\sum_{i=1}^{n}[\sigma+H_{i}(\frac{X_{i}\omega-y_{i}}{\sigma}) \sigma]+\alpha||\omega||_{2}^{2} \tag{13}\]  
Figure 2: Variation of BRDF in the red (**a,b**) and infrared bands (**c,d**) in PP and CPP planes with the observation zenith angle. \[H_{i}= \begin{cases}\dfrac{1}{2}[\rho_{i}-f(\rho_{i}^{\prime})]^{2},&|\rho_{i}-f( \rho_{i}^{\prime})|\leq\delta\\ \delta|\rho_{i}-f(\rho_{i}^{\prime})|-\dfrac{1}{2}\delta^{2},&|\rho_{i}-f(\rho_ {i}^{\prime})|>\delta\end{cases}  \tag{14}\]  
The parameter \(X\) is the training vector, which is the simulated land surface bidirectional reflectance obtained by the BRDF archetype under the same geometric conditions as the actual observation, in the form of an n-dimensional one-element array; \(y\) is the target vector, which is the land surface bidirectional reflectance of the land surface observed by the satellite, in the form of a one-dimensional n-element array. \(\omega\) is the coefficient vector, \(\alpha\) is the complexity parameter, and \(\sigma^{2}\) is defined as the variance estimate of the noise. \(f(\rho_{i}^{\prime})\) is a linear function of \(\rho_{i}^{\prime}\). For the detailed algorithm solution of Huber regression in the above formula, see the book [32]. To obtain the optimal archetype coefficients and improve the regression fitting accuracy, define the fitting parameter epsilon as \(\left(\rho_{i}-f(\rho_{i}^{\prime})\right)/\delta\) and set the epsilon range to [1.0.1,1.1,2.1,35,1.5,1.75,2.0,2.5,3.0]. Since the Huber loss fitting regression involves many parameters, to make the reading clear, the definitions of the parameters are listed in Table 1 in detail. The sum of the regression loss of each BRDF archetype in different epsilon parameters is compared, and the adjustment coefficients A and B corresponding to the smallest loss are taken as the final BRDF archetype coefficients.