Retrieval of Leaf Area Index by Linking the PROSAIL and Ross-Li BRDF Models Using MODIS BRDF Data

: Canopy structure parameters (e.g., leaf area index (LAI)) are key variables of most climate and ecology models. Currently, satellite-observed reﬂectances at a few viewing angles are often directly used for vegetation structure parameter retrieval; therefore, the information content of multi-angular observations that are sensitive to canopy structure in theory cannot be sufﬁciently considered. In this study, we proposed a novel method to retrieve LAI based on modelled multi-angular reﬂectances at sufﬁcient sun-viewing geometries, by linking the PROSAIL model with a kernel-driven Ross-Li bi-directional reﬂectance function (BRDF) model using the MODIS BRDF parameter product. First, BRDF sensitivity to the PROSAIL input parameters was investigated to reduce the insensitive parameters. Then, MODIS BRDF parameters were used to model sufﬁcient multi-angular reﬂectances. By comparing these reference MODIS reﬂectances with simulated PROSAIL reﬂectances within the range of the sensitive input parameters in the same geometries, the optimal vegetation parameters were determined by searching the minimum discrepancies between them. In addition, a signiﬁcantly linear relationship between the average leaf angle (ALA) and the coefﬁcient of the volumetric scattering kernel of the Ross-Li model in the near-infrared band was built, which can narrow the search scope of the ALA and accelerate the retrieval. In the validation, the proposed method attains a higher consistency (root mean square error (RMSE) = 1.13, bias = − 0.19, and relative RMSE (RRMSE) = 36.8%) with ﬁeld-measured LAIs and 30-m LAI maps for crops than that obtained with the MODIS LAI product. The results indicate the vegetation inversion potential of sufﬁcient multi-angular data and the ALA relationship, and this method presents promise for large-scale LAI estimation.


Introduction
Vegetation monitoring is critical for the development of adaptation strategies to address the challenges caused by climate change and human activities in ecosystems, such as global warming and precision farming [1]. In particular, satellite remote sensing has provided an effective way to perform quick global investigations [2], where the vegetation structure and biochemical parameters commonly considered in monitoring are derived from physical canopy reflectance models. Among the vegetation parameters, the leaf area satellite observations at a few angles [7,8]. To utilize sufficient anisotropy reflectances, multiangle data with a good sampling distribution at 397 sun-viewing angles were modelled by the well-known MODIS BRDF parameter product and a hotspot-revised Ross-Li model [39,40]. Despite of the different limitations behind the model mechanism, a good BRDF agreement has been represented in linking the Ross-Li model and the PROSAIL canopy BRDF model [42], which includes several essential vegetation parameters [43,44]. Therefore, the robust PROSAIL model was selected for LAI retrieval from these modelled MODIS multi-angular data. First, the sensitivities of the BRDF variables to vegetation parameters of the PROSAIL model were investigated to reduce the free variables in the physical model inversion. Then, the estimated LAIs were validated based on field measurements and high-resolution LAI maps, and the reasons for the estimation error were analysed. Finally, some explanations, limitations, and future expectations were discussed.

Canopy BRDF Simulations with the PROSAIL Model
The robust canopy radiative transfer model, i.e., the PROSAIL model, was adopted for vegetation parameter inversion, considering its simple input parameters and good BRDF consistency with the Ross-Li model [42]. The latest 4SAIL canopy BRDF model considering the hotspot effect and the PROSPECT-5 model, which describes the leaf optical reflectance and transmittance, are referenced in this study (http://teledetection.ipgp. jussieu.fr/prosail/, accessed on 10 August 2021) [43][44][45]. This 1-D model adopts a simplified assumption that considers horizontal homogeneity, where N layers of leaves are contained in the model scene. The input model parameters are listed in Table 1, covering the main structure and biology information of the vegetation. First, the leaf reflectance ρ l and transmittance τ l can be calculated from leaf parameters, as given in Equation (1), which are then input into the SAIL model to calculate the canopy reflectance ρ c , as expressed in Equation (2). Following the general process, the PROSAIL model can simulate canopy reflectance spectra in the optical domain from 400-2500 nm at a 1-nm resolution in addition to BRFs in an arbitrary sun-viewing geometry, which have been widely applied in vegetation inversion [11,[43][44][45][46][47][48][49][50].
ρ l (λ) & τ l (λ) = PROSPECT(N s , C ab , C ar , C brown , C w , C m ) (1) ρ c (λ) = 4SAIL(ρ l (λ), τ l (λ), LAI, ALA, H spot , P soil , SKYL, SZA, VZA, ϕ) To support vegetation inversion, a comprehensive simulation dataset containing 20,000 vegetation parameter combinations was generated with the Saltelli periodic function [51] via uniform sampling (Figure 1a) of seven leaf and canopy parameters. The parameters varied within the ranges listed in Table 1, which references the experimental designs of previous studies [39,52,53] except for C ar , C brown , H spot , and SKYL. Considering that most leaves do not contain brown pigment [52] and that the intrinsic BRDF is not necessarily related to the diffuse irradiance component [54], these two parameters were set to zero. In addition, the insensitive parameters of C ar and H spot have typically been defined as constants in vegetation inversion [11,14,52]. Compared to a previous study [42], C ar and H spot were defined as constants, and more reasonable ranges were considered for C ab , LAI, and average leaf angle (ALA) based on the application studies mentioned above.
Subsequently, we simulated the corresponding directional reflectance values with the PROSAIL model by inputting the 20,000 sets of vegetation parameters. Most LAI products are retrieved from data in the red and NIR bands due to their high sensitivity [55]. Therefore, PROSAIL simulations were performed at the central wavelength in these two MODIS bands (645 and 858 nm). Moreover, 397 sun-viewing angles were selected to simulate the BRFs, as indicated in Table 1 and Figure 1b, which can suitably reflect the general BRDF variance [42,56]. The PROSAIL simulation dataset of the BRFs was applied in this study, including for the sensitivity analysis and LAI estimation.

Parameters Unit Common Value Search Range (Step)
Leaf scale leaf structure parameter (N s ) -1.

500-m MODIS BRDF Parameter Product
The 500-m MODIS global BRDF parameter products provide the weight coefficients of three typical scattering components in the Ross-Li model, including isotropic scattering, volumetric scattering, and geometric-optical scattering [20,21]. The RossThick-LiSpars-eReciprocal (RTLSR) model was chosen as the operational algorithm to generate the MODIS BRDF parameter product [13,18], and its general expression is given in Equation (3). Numerous validation studies have demonstrated the good accuracy of the Ross-Li model [57][58][59][60]. In Equation (3), R denotes the surface directional reflectance as a function of three angles (solar zenith angle (SZA) (θs), viewing zenith angle (VZA) (θv), and relative azimuth angle (RAA) (φ)) and the waveband, which comprises three scattering kernel types: isotropic scattering kernel (1.0), volumetric scattering kernel (Kvol), and geometricoptical scattering kernel (Kgeo). By inputting multi-angular reflectance data into the model, the three weight coefficients of fiso, fvol, and fgeo can be determined based on least-squares regression. Notably, the three coefficients were restricted to be no less than 0.0 to represent a positive contribution of different kinds of scattering. Subsequently, the directional reflectance at arbitrary orientations can be simulated.

500-m MODIS BRDF Parameter Product
The 500-m MODIS global BRDF parameter products provide the weight coefficients of three typical scattering components in the Ross-Li model, including isotropic scattering, volumetric scattering, and geometric-optical scattering [20,21]. The RossThick-LiSparseReciprocal (RTLSR) model was chosen as the operational algorithm to generate the MODIS BRDF parameter product [13,18], and its general expression is given in Equation (3). Numerous validation studies have demonstrated the good accuracy of the Ross-Li model [57][58][59][60]. In Equation (3), R denotes the surface directional reflectance as a function of three angles (solar zenith angle (SZA) (θ s ), viewing zenith angle (VZA) (θ v ), and relative azimuth angle (RAA) (ϕ)) and the waveband, which comprises three scattering kernel types: isotropic scattering kernel (1.0), volumetric scattering kernel (K vol ), and geometricoptical scattering kernel (K geo ). By inputting multi-angular reflectance data into the model, the three weight coefficients of f iso , f vol , and f geo can be determined based on least-squares regression. Notably, the three coefficients were restricted to be no less than 0.0 to represent a positive contribution of different kinds of scattering. Subsequently, the directional reflectance at arbitrary orientations can be simulated.

of 24
The latest Collection V006 daily BRDF parameter product (i.e., MCD43A1) and quality data (i.e., MCD43A2) [20] (https://urs.earthdata.nasa.gov/, accessed on 10 August 2021) were implemented in this study, which are concurrent with the field LAI measurements and LAI maps introduced in Section 2.3. MCD43A1 data were retrieved from the RTLSR model by inputting 16-day aggregation values of the multi-angular reflectance of the MODIS sensors onboard the Terra (overpass occurs at 10:30 am) and Aqua (overpass occurs at 1:30 pm) satellites. BRDF parameters were retrieved with the main algorithm (full inversion, quality flag = 0/1) and backup algorithm (magnitude inversion, quality flag = 2/3) in the red and NIR bands, respectively. High-quality BRDF data were retrieved via the Ross-Li model with good angular sampling and small fitting root mean square errors (RMSEs) in the inversion results. When a high-quality full inversion could not be accomplished due to poor angular sampling or insufficient input observations (i.e., obs. num. < 7), a priori knowledge method of the backup algorithm was employed (also called a magnitude inversion), which performs quite well in many situations [20,21,59,61,62]. The backup algorithm relies on an a priori archetypal BRDF shape with each pixel on the globe, and then the satellite-observed directional reflectances are used to adjust these archetypal shapes by a multiplicative factor to obtain suitable BRDF variances at each pixel [20,21].
These MODIS BRDF parameters were applied to model BRFs with angle sampling, as shown in Figure 1b, using the hotspot-revised Ross-Li model [39,40], and these modelled MODIS BRFs then functioned as reference data to search for the best vegetation estimation among the 20,000 PROSAIL vegetation parameter datasets, as introduced in Section 2.1.

LAI Measurements at the 500-m Plot Level and with 30-m LAI Maps
Field-measured LAI values at the 500-m plot level and 30-m LAI maps were adopted for method validation, which encompassed two sites (Honghe and Hailun) in Heilongjiang Province, northeastern China.
There were 180 sets in total of field-measured LAI data, which were collected by Fang et al. in 2012Fang et al. in , 2013, and 2016 (https://doi.pangaea.de/10.1594/PANGAEA.900090, accessed on 10 August 2021) [55,63,64]. Both sites cover a homogeneous area of approximately 30 km 2 , and 5 plots were selected to perform measurements at each site. The first site is located at Honghe Farm (47 •  Field measurements were carried out under twilight and overcast conditions to avoid direct sun flecks, from shortly after the crop greening stage to the harvest readiness stage. LAI-2200 measurements were used in this study, which agreed well with the results obtained via the destructive leaf sampling method in most of the growing seasons from DOYs 160 to 230. These LAI measurements contained the main crop types, and they have been considered in the validation of several essential satellite LAI products, such as the MODIS LAI product generated by the main algorithm [55]. These field-measured LAI values were used to validate the retrieval accuracy of the new method proposed in this study, as described in Sections 3.3, 4.2 and 4.3.
In addition, LAI maps concurrent with these LAI measurements with a spatial resolution of 30 m were applied for method validation, which were provided by a previous study [55]. High-resolution reference LAI maps were derived from HJ-1, Landsat 7, and Sentinel-2A images, which agree well with the field-measured LAI values. Moreover, these maps were adopted in the validation of the moderate-resolution LAI products mentioned above, where average values of 3 × 3 pixels were considered to cross validate the moderateresolution LAI products and concurrent reference LAI maps. In this study, 284 sets of aggregated LAI values at 1.5 km derived from 30-m LAI maps were used to cross validate our proposed novel method, as introduced in Section 4.4, and the data numbers were 108 and 176 in Honghe and Hailun, respectively.

Methods
Considering the satellite limitations when capturing BRDF data in the viewing hemisphere as well as their underlying observation noise, the semi-empirical Ross-Li BRDF model has been widely adopted to model BRDF data at arbitrary sun-viewing geometries of good quality from satellite-observed reflectances at a few viewing angles. In this study, a novel method for the simultaneous retrieval of the LAI and other vegetation parameters based on sufficient multi-angular reflectance data using the Ross-Li-derived MODIS BRDF parameter product was developed, in conjunction with the hotspot-revised Ross-Li model and physical PROSAIL BRDF model, considering their good BRDF consistency [42]. The new method mainly includes two parts: sensitivity analysis and vegetation parameter retrieval, and the flowchart is shown as Figure 2.
First, BRDF sensitivity to the PROSAIL input parameters was investigated to reduce the number of insensitive parameters in the vegetation inversion. The analysis was performed based on PROSAIL-simulated multi-angular reflectance data at 397 angles in the red and NIR bands, as described in Section 2.1, as well as the corresponding BRDF variables retrieved with the Ross-Li model, including the three kernel coefficients and anisotropy flat index (AFX) [65]. To maintain consistency with the MODIS BRDF parameter product used for LAI retrieval in this study, the RTLSR model was employed. Subsequently, the extended Fourier amplitude sensitivity test (EFAST) method was implemented to explore the sensitivity of the BRDF variance to the vegetation biophysical/structure parameters, thereby further investigating the empirical relationship between these highly sensitive parameters. Encouragingly, a significantly linear relation between the ALA and f vol in the NIR band was established.

Sensitivity Analysis of the BRDF to Vegetation Parameters
Previous studies have indicated that compensation occurs between canopy variables such as the LAI, leaf inclination, Cab, and leaf water content [45,66,67], and these mutually dependent parameters can therefore be hard to accurately determine during retrieval. Via the application of the obtained sensitivity analysis results, simultaneous retrieval was accomplished to estimate the sensitive vegetation parameters, including the LAI, from the MODIS BRDF parameter product by linking the PROSAIL and Ross-Li BRDF models. BRFs with the same 397 angles of the PROSAIL simulations, as shown in Figure 1b, were modelled using the MODIS BRDF parameters and Ross-Li model, and the latest hotspotcorrected version of the Ross-Li model was employed here to better simulate the hotspot reflectance [39,40]. Then, these MODIS multi-angular reflectance values functioned as reference data to retrieve the optimal vegetation parameters by simultaneously searching the minimum discrepancies and the PROSAIL-simulated reflectances from the PROSAIL vegetation parameter database. Both wide and local searches with the empirically calculated ALA were designed. Finally, the reference LAI data in Section 2.3 were utilized for method validation.

Sensitivity Analysis of the BRDF to Vegetation Parameters
Previous studies have indicated that compensation occurs between canopy variables such as the LAI, leaf inclination, C ab , and leaf water content [45,66,67], and these mutually dependent parameters can therefore be hard to accurately determine during retrieval. Therefore, the sensitivity of the BRDF variance to vegetation parameters was investigated to reduce the number of insensitive and dependent parameters. The vegetation parameters refer to the simulated 20,000 PROSAIL parameter database, and the BRDF variables were retrieved via the Ross-Li model corresponding to the simulated PROSAIL BRFs. Considering that the MODIS BRDF parameter product used for LAI retrieval was generated from the RTLSR model, we further applied this model to calculate the BRDF variables in the sensitivity analysis.
In detail, the directional reflectance values obtained from the PROSAIL model in Section 2.1 were then input into the RTLSR model to invert the three kernel coefficients (f iso , f vol , and f geo ) and BRDF index (AFX). The AFX was calculated as the ratio of the white sky albedo (WSA) to f iso [65], which is sensitive to the scattering type. An AFX higher than 1.0 indicates that volumetric scattering is prominent, while an AFX lower than 1.0 indicates dominant geometric-optical scattering. Based on the 20,000 sets of BRDF variables and concurrent vegetation parameters, their sensitivity was investigated via the EFAST statistic method and linear regression analysis.
The EFAST method was first used for the sensitivity analysis [68], which provides the global sensitivity to evaluate the contribution of the input parameters to the output by considering the interactions among the input data. The calculation process is given as follows: first, the total variance in the model output attributed to the input parameters was divided into two parts: the variance resulting from the specific X i (D i ) value and the variance resulting from the interactions between the input parameters under multiple combinations. Then, the local sensitivity S i and mutual sensitivity S i . . . k could be obtained by dividing the total variance D by the corresponding individual and interaction variances, respectively. Finally, the global sensitivity S i T could be obtained as the accumulation of the local sensitivity S i and all the mutual sensitivity values related to X i . In total, the global sensitivity of the four BRDF variables (f iso , f vol , f geo , and AFX) to seven vegetation parameters (N s , C ab , C w , C m , LAI, ALA, and P soil ) was calculated to perform a qualitative analysis.
Remote Sens. 2021, 13, 4911 8 of 24 To further explore the quantitative relations between the BRDF variables and vegetation parameters, the highly sensitive parameters determined according to the EFAST method were used for linear regression. Notably, only parameters that met the above requirement were considered in the sensitivity analysis, where the fitting RMSE values between the PROSAIL BRFs and BRFs inverted from the RTLSR model were smaller than 0.02 and 0.05 in the red and NIR bands, respectively [42]. To reduce the influence of data uncertainty, the 10-fold cross-validation method was utilized in the regression, and 90% of the data was selected for regression, and the remaining 10% was reserved for validation. In total, 10 independent regressions were performed, and those linear relationships with the minimum regression error were recorded as the final result. Encouragingly, a significantly linear relationship between the ALA and f vol in the NIR band was established, which was used to develop the LAI retrieval method.

LAI Estimation from MODIS BRDF Data by Linking the PROSAIL and Ross-Li Models
Based on the sensitivity analysis of the BRDF sensitivity to vegetation parameters in Section 3.1, a general vegetation inversion method was proposed in this study to retrieve these sensitive vegetation parameters. The MODIS BRDF parameter product was applied for multi-angular reflectance modelling at sufficient sun-viewing geometries with the Ross-Li model, which was then used to retrieve the vegetation parameters from the robust PROSAIL physical BRDF model.
First, the modelled BRFs from the MODIS BRDF parameters were employed as reference data in the vegetation parameter inversion. BRFs considering 397 angles, as shown in Figure 1b, were modelled with the Ross-Li model based on the daily MODIS BRDF parameters introduced in Section 2.2, and a hotspot-revised version of the Ross-Li model, namely, RTLSR_C, was used here [39,40]. The revised model achieves a higher hotspot capability than that of the RTLSR model by applying the corrected exponential hotspot function developed by Chen and Cihlar [69] to the volumetric and geometric-optical scattering kernels. The best hotspot parameters of the height (C 1 ) and width (C 2 ) of the RTLSR_C model in several typical bands were given based on a wide search of sufficient hotspot data, as listed in Table 2 [39], and the hotspot parameters at 645 and 858 nm for MODIS (C 1 = 0.5, C 2 = 3.4 • at 645 nm; C 1 = 0.5, C 2 = 3.0 • at 858 nm) were used to model the MODIS BRFs based on the obtained MODIS BRDF parameters.
Then, two search strategies, including wide and local search strategies based on the f vol -ALA relationship, were developed to retrieve the optimal vegetation parameters from the 20,000 vegetation parameter database introduced in Section 2.1. The cost function is expressed below and calculates the discrepancies between the simulated PROSAIL BRFs and reference MODIS BRFs, where ρ i reference and ρ i PROSAIL denote the simulated MODIS BRFs and PROSAIL BRFs, respectively, for a specific group of PROSAIL input parameters among the 20,000 vegetation parameter databases. Considering the BRFs in both the red and NIR bands, the value of N is 397 times two (i.e., 794). Finally, the first 50 records with the minimum cost values were selected, and the corresponding LAI values were averaged as the retrieved LAI result.
In the wide search strategy, the optimal estimation is determined by traversing all 20,000 vegetation parameter datasets. In the local search strategy with the ALA, the linear relationship between the ALA and f vol in the NIR band mentioned in Section 3.1 was applied, and the empirical ALA could thus be obtained from the MODIS f vol value in the NIR band. Consequently, this narrowed the search scope within the theoretical error range of the empirical ALA rather than the whole reasonable range, which could greatly improve the search speed.

LAI Validation and Error Analysis
The 180 LAI measurements at the 500-m plot level and the 30-m LAI maps covering the main crop types were utilized for method validation, as introduced in Section 2.3. The evaluation indicators include the RMSE, relative RMSE (RRMSE), bias, and coefficient of determination (R 2 ) between the reference LAIs and retrieved data, where the RMSE and RRMSE are expressed below. LAI i obs and LAI i estimate denote the reference LAI and the estimated LAI with the PROSAIL model, respectively.
Then, we investigated the error underlying the proposed method to further improve the vegetation estimation accuracy. By comparing the LAI estimation accuracies based on the MODIS BRDF data between the main and backup algorithms, we investigated the influence of the MODIS BRDF quality on the LAI estimation accuracy. As expected, BRDF data from the main algorithm have better LAI estimations than those from the backup algorithm. Therefore, the BRDF data retrieved with the backup algorithm were alternated by those retrieved with the main algorithm at the nearest time to reduce the impact of the BRDF quality on the proposed method. In addition, the LAI estimation accuracy of the new method proposed in this study was cross compared to published accuracy results for MODIS LAI products generated from satellite-observed reflectances at a few viewing angles [55]. Considering the effects of the point spread function and geometric distortion of satellite data, 30-m LAI maps were further used for cross validation in the upscaling process to obtain a 1.5 km × 1.5 km area [55].

Analysis Based on EFAST
Based on the 20,000 sets of PROSAIL multi-angular canopy reflectance simulations and BRDF variables retrieved from the Ross-Li RTLSR model, the global sensitivity of the BRDF to the vegetation parameters was calculated, as shown in Figure 3. The sensitivity of the nadir reflectance to the vegetation parameters has been investigated in previous studies [52]. This study first analysed the results regarding BRDF indicators. Figure 3 shows that most BRDF variables are highly sensitive to many vegetation parameters in the red band (N s , C ab , LAI, ALA and Psoil). In the NIR band, N s , LAI, and ALA appear to be the most sensitive parameters, and f vol reveals the most significant sensitivity to the ALA among all the statistics. It is foreseeable for the high sensitivity of canopy BRDF to the two essential canopy structure parameters of LAI and ALA, which have shown similar results associated to nadir reflectance [45,52]. In addition, the sensitivity to C ab may be caused by high chlorophyll absorption in the red band, and the non-negligible BRDF change in the background soil associated with Psoil also makes important contributions to the overall canopy BRDF [70]. The leaf structure parameter N s represents the number of compact layers specifying the average number of air/cell wall interfaces within the mesophyll and the leaf biochemical content, which changes with vegetation species and phenology that are sensitive to canopy BRDF [42,45]. N s for the monocotyledons changes from 1.0 to 1.5, while values for the dicotyledons range from 1.5 to 2.5. Meanwhile, the senescent leaves always show an N s larger than 2.5.
These results generally comply with the sensitivity of the canopy nadir reflectance to PROSAIL variables [45,52], although there are few connections between canopy reflectance anisotropy and some vegetation parameters, especially for the chlorophyll and water content in the NIR band. Nevertheless, all seven vegetation parameters have showed evident sensitivities in the red or NIR bands, which can be considered to be an overall sensitivity to the seven vegetation parameters. In addition, observations in the two bands have been widely utilized together to generate a series of LAI products [55], and thus all seven PROSAIL variables were applied in the vegetation inversion algorithm owing to their sensitivity.
ALA among all the statistics. It is foreseeable for the high sensitivity of canopy BRDF to the two essential canopy structure parameters of LAI and ALA, which have shown similar results associated to nadir reflectance [45,52]. In addition, the sensitivity to Cab may be caused by high chlorophyll absorption in the red band, and the non-negligible BRDF change in the background soil associated with Psoil also makes important contributions to the overall canopy BRDF [70]. The leaf structure parameter Ns represents the number of compact layers specifying the average number of air/cell wall interfaces within the mesophyll and the leaf biochemical content, which changes with vegetation species and phenology that are sensitive to canopy BRDF [42,45]. Ns for the monocotyledons changes from 1.0 to 1.5, while values for the dicotyledons range from 1.5 to 2.5. Meanwhile, the senescent leaves always show an Ns larger than 2.5.
These results generally comply with the sensitivity of the canopy nadir reflectance to PROSAIL variables [45,52], although there are few connections between canopy reflectance anisotropy and some vegetation parameters, especially for the chlorophyll and water content in the NIR band. Nevertheless, all seven vegetation parameters have showed evident sensitivities in the red or NIR bands, which can be considered to be an overall sensitivity to the seven vegetation parameters. In addition, observations in the two bands have been widely utilized together to generate a series of LAI products [55], and thus all seven PROSAIL variables were applied in the vegetation inversion algorithm owing to their sensitivity.

Modelling and Evaluation of Linear Relationships
The sensitive parameters determined with the EFAST method were further utilized for modelling and evaluation of their linear relations, where only the data meeting the fit-RMSE requirements of the Ross-Li model in the red and NIR bands were selected. Among the 20,000 sets of PROSAIL canopy BRFs and corresponding Ross-Li BRDF variables, there are 15,707 sets of data where the fit-RMSE is less than 0.02 and 0.05 in the red and NIR bands, respectively. When performing 10-fold cross regression and validation, all 15,707 datasets were divided into 10 parts through equal-interval sampling.
Encouragingly, a significantly linear relationship between the ALA and fvol in the NIR band was obtained. The regression coefficients and evaluation accuracies considering the 10 parts are listed in Table 3. The table indicates that the fit coefficients are extremely similar for all 10 data groups, exhibiting a highly significant F at the significance level of 0.01 in the F test. The equations of the F test are shown in (11)-(13), where F is calculated as the ratio of the averaged regression sum of squares (SSR) to the averaged error sum of squares (SSE). The variables yi, ӯ, and ŷi refer to the values of the ith ALA, the averaged ALAs, and the predicted ALAs by the linear regression relation, respectively. Considering that only a small amount of data is sparsely distributed in the ALA evaluation, the accuracy indicators were calculated according to all the evaluation data as well as the highdensity data, except the sparse evaluation data. By considering the errors in both regression and evaluation, the fourth regression showed a relatively smaller error for both all

Modelling and Evaluation of Linear Relationships
The sensitive parameters determined with the EFAST method were further utilized for modelling and evaluation of their linear relations, where only the data meeting the fit-RMSE requirements of the Ross-Li model in the red and NIR bands were selected. Among the 20,000 sets of PROSAIL canopy BRFs and corresponding Ross-Li BRDF variables, there are 15,707 sets of data where the fit-RMSE is less than 0.02 and 0.05 in the red and NIR bands, respectively. When performing 10-fold cross regression and validation, all 15,707 datasets were divided into 10 parts through equal-interval sampling.
Encouragingly, a significantly linear relationship between the ALA and f vol in the NIR band was obtained. The regression coefficients and evaluation accuracies considering the 10 parts are listed in Table 3. The table indicates that the fit coefficients are extremely similar for all 10 data groups, exhibiting a highly significant F at the significance level of 0.01 in the F test. The equations of the F test are shown in (11)- (13), where F is calculated as the ratio of the averaged regression sum of squares (SSR) to the averaged error sum of squares (SSE). The variables y i , y , andŷ i refer to the values of the ith ALA, the averaged ALAs, and the predicted ALAs by the linear regression relation, respectively. Considering that only a small amount of data is sparsely distributed in the ALA evaluation, the accuracy indicators were calculated according to all the evaluation data as well as the high-density data, except the sparse evaluation data. By considering the errors in both regression and evaluation, the fourth regression showed a relatively smaller error for both all data and high-density data. Therefore, the results from the fourth data group were recorded as the final regression results. The final optimal linear relation is shown in Equation (14), and the fitting results are shown in the density scatter plots (Figure 4). The scatter plots are sliced into five equaldensity levels from 1 (highest) to 5 (lowest), by classifying the data number with intervals of 1.0 • and 0.01 for ALA and f vol , respectively. Considering the effective range of the ALA from 10 • to 85 • , as explained in Table 1, only f vol in the NIR band smaller than 0.3813 could be used to retrieve the empirical ALA based on this linear relationship. Figure 4 clearly shows that only a small amount of data is sparsely distributed in both the regression and ALA evaluation, and most high-density data are clustered. The high R 2 value of 0.98 and significance F value, as shown in Figure 4a, indicate that the linear relationship achieves a high quality. Similarly, independent evaluation also reveals a high R 2 value between the reference ALA and estimated ALA based on 10% of the 15,707 datasets, which is 0.87 for the total data and 0.99 for the high-density data. By removing a small amount of sparse data, a generally good accuracy was attained for empirical ALA retrieval based on Equation (14) (RMSE = 2.07 • , bias = −0.09 • , R 2 = 0.99). The significant linear relationship may be attributed to the intrinsic radiative transfer hypotheses of the PROSAIL model, where multiple scattering (i.e., volumetric scattering) can easily arise in horizontal homogeneous scenes [45]. In addition, the relationship in the NIR band complies with a previous study, which indicates that volumetric scattering is prominent in the NIR band for the PROSAIL model [42].

Validation of the LAI Estimations Based on Field Measurements
A total of 180 LAI measurements at the 500-m plot level were used for algorithm validation, including wide and local search strategies via the empirical ALA based on Equation (14). Among the 180 sets of field-measured LAI values, nine sets of concurrent MODIS BRDF parameters comprised fill values, for which BRDF data retrieved with the main and backup algorithms at the nearest time were added as alternative data to guarantee data integrity. There are 103 sets of BRDF data that yielded reasonable ALA values ranging from 10° to 85°, while the remaining data could not employ the local search method because the fvol value in the NIR band was larger than 0.3813. Considering that the theoretical estimation-based RMSE of the linear estimated ALA is 2.07°, vegetation parameter combinations with ALAs varying by ±3° around the empirical value were used in the local search strategy. Notably, the local search range directly decreased to approximately 1600 sets of data, while all 20,000 sets of data were traversed in the wide search. As a result, the local search was 11.5 times faster than the wide search, which greatly improved the retrieval speed. Figure 5 and Table 4 illustrate the validation results of the estimated LAIs using the proposed method based on LAI field measurements. During the late growing season after DOY 240 (bluish points), the estimated LAIs were lower than the reference LAIs, which agrees with the validation results of the MODIS LAI product [55]. The local search method based on the empirical ALA is not applicable mainly to paddy rice in Honghe, where fvol in the NIR band for 70 sets of data was larger than 0.3813 among all 110 sets. Regarding the crops in Hailun, only in 7 sets of data did the retrieved empirical ALA occur beyond the range in the 70 sets. With the use of the same data, the LAI estimation accuracy was obviously improved with the local search method (the third column of Figure 5) over the wide search method (the second column of Figure 5). Consequently, fusion results could be obtained by combining all of the results obtained in the local search and the remaining results of the wide search. As expected, the fused LAI retrievals of these two methods (the fourth column of Figure 5; RMSE = 1.34, bias = −0.52, R 2 = 0.49, RRMSE = 41.3%) exhibited a higher LAI estimation accuracy than that achieved with the wide search method (the first column of Figure 5), as well as the MODIS LAI product generated with the main algorithm (RMSE = 1.57, bias = −0.21, R 2 = 0.24, RRMSE = 48.3%) [8] given the same data, as listed in Table 4. The results highlight the importance of sufficient BRDF information in vegetation parameter retrieval over the utilization of satellite observation data at a few viewing angles, such as the MODIS LAI product.

Validation of the LAI Estimations Based on Field Measurements
A total of 180 LAI measurements at the 500-m plot level were used for algorithm validation, including wide and local search strategies via the empirical ALA based on Equation (14). Among the 180 sets of field-measured LAI values, nine sets of concurrent MODIS BRDF parameters comprised fill values, for which BRDF data retrieved with the main and backup algorithms at the nearest time were added as alternative data to guarantee data integrity. There are 103 sets of BRDF data that yielded reasonable ALA values ranging from 10 • to 85 • , while the remaining data could not employ the local search method because the f vol value in the NIR band was larger than 0.3813. Considering that the theoretical estimation-based RMSE of the linear estimated ALA is 2.07 • , vegetation parameter combinations with ALAs varying by ±3 • around the empirical value were used in the local search strategy. Notably, the local search range directly decreased to approximately 1600 sets of data, while all 20,000 sets of data were traversed in the wide search. As a result, the local search was 11.5 times faster than the wide search, which greatly improved the retrieval speed. Figure 5 and Table 4 illustrate the validation results of the estimated LAIs using the proposed method based on LAI field measurements. During the late growing season after DOY 240 (bluish points), the estimated LAIs were lower than the reference LAIs, which agrees with the validation results of the MODIS LAI product [55]. The local search method based on the empirical ALA is not applicable mainly to paddy rice in Honghe, where f vol in the NIR band for 70 sets of data was larger than 0.3813 among all 110 sets. Regarding the crops in Hailun, only in 7 sets of data did the retrieved empirical ALA occur beyond the range in the 70 sets. With the use of the same data, the LAI estimation accuracy was obviously improved with the local search method (the third column of Figure 5) over the wide search method (the second column of Figure 5). Consequently, fusion results could be obtained by combining all of the results obtained in the local search and the remaining results of the wide search. As expected, the fused LAI retrievals of these two methods (the fourth column of Figure 5; RMSE = 1.34, bias = −0.52, R 2 = 0.49, RRMSE = 41.3%) exhibited a higher LAI estimation accuracy than that achieved with the wide search method (the first column of Figure 5), as well as the MODIS LAI product generated with the main algorithm (RMSE = 1.57, bias = −0.21, R 2 = 0.24, RRMSE = 48.3%) [8] given the same data, as listed in Table 4. The results highlight the importance of sufficient BRDF information in vegetation parameter retrieval over the utilization of satellite observation data at a few viewing angles, such as the MODIS LAI product. Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 23  There are five methods included in total: (1) MODIS refers to the latest Collection V006 MODIS LAI product (i.e., MCD15A2H), where only data generated by the main  There are five methods included in total: (1) MODIS refers to the latest Collection V006 MODIS LAI product (i.e., MCD15A2H), where only data generated by the main algorithm are used [55]; (2) wide search refers to results retrieved by the wide research method; (3) wide search of local search data refers to results retrieved with the wide research method considering the same LAI data as those considered under the local search method; (4) local search by the ALA refers to results retrieved with the local research method; (5) fused results refers to the combined results obtained from the local search method and the remaining results of the wide search method.

Error Analysis and Method Improvement
We further analyzed the reason for the estimation error, especially in the quality of the MODIS BRDF parameters. The ALA validation results for the MODIS BRDF parameters of different qualities are shown in Figure 6, and the estimation accuracy based on high-quality BRDF data (QA = 0/1) is obviously higher than that based on data retrieved from the backup algorithm (QA = 2/3). The difference in LAI estimation accuracy resulting from the BRDF quality is easy to understand because the backup algorithm is employed when the main algorithm cannot be implemented due to poor angular sampling or insufficient input observations [20,21,59,61,62]. The analysis further emphasizes the significance of sufficient high-quality BRDF information for vegetation parameter retrieval.

Error Analysis and Method Improvement
We further analyzed the reason for the estimation error, especially in the quality of the MODIS BRDF parameters. The ALA validation results for the MODIS BRDF parameters of different qualities are shown in Figure 6, and the estimation accuracy based on high-quality BRDF data (QA = 0/1) is obviously higher than that based on data retrieved from the backup algorithm (QA = 2/3). The difference in LAI estimation accuracy resulting from the BRDF quality is easy to understand because the backup algorithm is employed when the main algorithm cannot be implemented due to poor angular sampling or insufficient input observations [20,21,59,61,62]. The analysis further emphasizes the significance of sufficient high-quality BRDF information for vegetation parameter retrieval. To reduce the impact of the BRDF quality on the proposed method and further improve the LAI retrieval accuracy, BRDF data obtained from the backup algorithm were tentatively alternated by data retrieved from the main algorithm at the nearest time. This strategy follows the assumption that there exists a similar BRDF pattern in half to all of the months as the MODIS [20] and POLDER BRDF products [22]. The updated validation results of the estimated LAI values using high-quality MODIS BRDF data are shown in Figure 7 and Table 5, which indicate a better accuracy than that using BRDF data with both the main and backup algorithms in Section 4.2. Similarly, the local search method via the ALA attains better estimations than does the wide search method for the same data, and the fused results yield the optimal LAI estimations (RMSE = 1.28, bias = −0.46, R 2 = 0.45, RRMSE = 39.5%) among all the methods proposed in this study. Unfortunately, estimation cannot be accomplished by the local search method considering more data than before. Similarly, the local search method still failed for mainly paddy rice in Honghe, but the data number increased from 70 to 93. For the crops in Hailun, the number of data points that failed to implement the local search method increased from 7 to 11. To reduce the impact of the BRDF quality on the proposed method and further improve the LAI retrieval accuracy, BRDF data obtained from the backup algorithm were tentatively alternated by data retrieved from the main algorithm at the nearest time. This strategy follows the assumption that there exists a similar BRDF pattern in half to all of the months as the MODIS [20] and POLDER BRDF products [22]. The updated validation results of the estimated LAI values using high-quality MODIS BRDF data are shown in Figure 7 and Table 5, which indicate a better accuracy than that using BRDF data with both the main and backup algorithms in Section 4.2. Similarly, the local search method via the ALA attains better estimations than does the wide search method for the same data, and the fused results yield the optimal LAI estimations (RMSE = 1.28, bias = −0.46, R 2 = 0.45, RRMSE = 39.5%) among all the methods proposed in this study. Unfortunately, estimation cannot be accomplished by the local search method considering more data than before. Similarly, the local search method still failed for mainly paddy rice in Honghe, but the data number increased from 70 to 93. For the crops in Hailun, the number of data points that failed to implement the local search method increased from 7 to 11. Remote Sens. 2021, 13, x FOR PEER REVIEW 15 of 23  The four methods are the same as those in Table 4.

General Flowchart for LAI Estimation and Validation Based on LAI Maps
According to the validation and error analysis results for the proposed method based on field measurements, a general flowchart for LAI retrieval by using the novel method  The four methods are the same as those in Table 4.

General Flowchart for LAI Estimation and Validation Based on LAI Maps
According to the validation and error analysis results for the proposed method based on field measurements, a general flowchart for LAI retrieval by using the novel method is shown in Figure 8. The process describes the flow scheme to estimate LAI via wide and local search strategies based on the MODIS BRDF parameter product.
Before method application, the data and methods should be prepared, including the PROSAIL model, hotspot-corrected Ross-Li RTLSR_C model, 20,000 PROSAIL vegetation parameter datasets, 397 angle samplings, and empirical function between the ALA and f vol in the NIR band, as expressed in Equation (14).
After preparation, we can start the estimation process. First, the quality of the MODIS BRDF parameters should be examined, and lower-quality BRDF data are alternated by high-quality data at the nearest time. Then, a second evaluation is performed to determine whether the value of f vol in the NIR band is less than 0.3813. If f vol satisfies this condition, the local search method is applied to estimate the final LAI, which searches the 20,000 vegetation parameter database only for those values approaching the empirical ALA calculated by the f vol -ALA relationship. In contrast, when the value of f vol exceeds the critical value, the wide search method is applied to estimate the LAI by searching the optimal results from the total 20,000 dataset. In this way, a continuous LAI estimation both in time and space can be obtained by jointly applying these two search methods. is shown in Figure 8. The process describes the flow scheme to estimate LAI via wide and local search strategies based on the MODIS BRDF parameter product. Before method application, the data and methods should be prepared, including the PROSAIL model, hotspot-corrected Ross-Li RTLSR_C model, 20,000 PROSAIL vegetation parameter datasets, 397 angle samplings, and empirical function between the ALA and fvol in the NIR band, as expressed in Equation (14).
After preparation, we can start the estimation process. First, the quality of the MODIS BRDF parameters should be examined, and lower-quality BRDF data are alternated by high-quality data at the nearest time. Then, a second evaluation is performed to determine whether the value of fvol in the NIR band is less than 0.3813. If fvol satisfies this condition, the local search method is applied to estimate the final LAI, which searches the 20,000 vegetation parameter database only for those values approaching the empirical ALA calculated by the fvol-ALA relationship. In contrast, when the value of fvol exceeds the critical value, the wide search method is applied to estimate the LAI by searching the optimal results from the total 20,000 dataset. In this way, a continuous LAI estimation both in time and space can be obtained by jointly applying these two search methods. Following the flowchart, 500-m LAI values estimated from the MODIS BRDF parameters concurrent with the high-resolution LAI maps were calculated. Comparisons of the retrieved LAI to the upscaled 30-m reference LAI maps at 500-m and 1.5-km resolutions (3 × 3 MODIS pixels) are shown in Figure 9. Generally, the estimated LAIs based on the proposed method present a variance trend similar to that of the reference data throughout the whole growing season, especially in Hailun, while the estimations are lower than the reference LAI data at the senescence stage, as shown in Figure 9a,c. Based on the scatter plots, an overall underestimation can be observed for paddy rice in Honghe, which agrees with the validation analysis at the 500-m plot level. For the maize, soybean, and sorghum crops in Hailun, a better fitting line close to the 1:1 line was obtained. Finally, the total Following the flowchart, 500-m LAI values estimated from the MODIS BRDF parameters concurrent with the high-resolution LAI maps were calculated. Comparisons of the retrieved LAI to the upscaled 30-m reference LAI maps at 500-m and 1.5-km resolutions (3 × 3 MODIS pixels) are shown in Figure 9. Generally, the estimated LAIs based on the proposed method present a variance trend similar to that of the reference data throughout the whole growing season, especially in Hailun, while the estimations are lower than the reference LAI data at the senescence stage, as shown in Figure 9a,c. Based on the scatter plots, an overall underestimation can be observed for paddy rice in Honghe, which agrees with the validation analysis at the 500-m plot level. For the maize, soybean, and sorghum crops in Hailun, a better fitting line close to the 1:1 line was obtained. Finally, the total accuracy for the crops at both sites (RMSE = 1.13, bias = −0.19, R 2 = 0.64, RRMSE = 36.8%) was higher than that at the plot level in Section 4.3 (RMSE = 1.28, bias = −0.46, R 2 = 0.45, RRMSE = 39.5%). Moreover, the estimated LAI values also indicate a better consistency with the reference data than does the MODIS LAI product (RMSE = 1.50, bias = −0.23, R 2 = 0.24, RRMSE = 47.1%) [55]. Similarly, there exists a better agreement with the LAI maps for the retrieved LAIs using the flowchart in this study than that for the MODIS LAI product, as described in Sections 4.2 and 4.3, which further indicates that modelled multi-angular reflectance data at sufficient sun-viewing geometries can provide effective BRDF information and thus can be applied in vegetation parameter retrieval.  [55]. Similarly, there exists a better agreement with the LAI maps for the retrieved LAIs using the flowchart in this study than that for the MODIS LAI product, as described in Sections 4.2 and 4.3, which further indicates that modelled multi-angular reflectance data at sufficient sun-viewing geometries can provide effective BRDF information and thus can be applied in vegetation parameter retrieval.

Discussion
The intrinsic surface anisotropic reflectance (i.e., the BRDF effect) is sensitive to vegetation structure parameters, and satellite-observed reflectance data can in turn promote global vegetation parameter retrieval. Nevertheless, BRDF information retrieved from satellite observations at a few viewing angles cannot sufficiently be considered in general. Therefore, modelled multi-angular reflectance values at sufficient sun-viewing geometries were used in this study, which were calculated based on the kernel-driven Ross-Li BRDF model by inputting the MODIS BRDF parameter product. As a result, compared to

Discussion
The intrinsic surface anisotropic reflectance (i.e., the BRDF effect) is sensitive to vegetation structure parameters, and satellite-observed reflectance data can in turn promote global vegetation parameter retrieval. Nevertheless, BRDF information retrieved from satellite observations at a few viewing angles cannot sufficiently be considered in general. Therefore, modelled multi-angular reflectance values at sufficient sun-viewing geometries were used in this study, which were calculated based on the kernel-driven Ross-Li BRDF model by inputting the MODIS BRDF parameter product. As a result, compared to direct satellite observations, these modelled, multi-angular data encompass sufficient angles, by which BRDF information can be better considered in theory. Moreover, cloud cover and other reasons often lead to noise in satellite observations, which may, in turn, affect the retrieval accuracy of vegetation parameters by directly using these observed data. However, the widely employed Ross-Li model for BRDF reconstruction (e.g., using the global MODIS BRDF parameter product) provides an advantage in terms of greatly smoothing random noise in satellite observations by theoretically considering the overall BRDF variance [13,18].
Based on these advantages, the potential of sufficient modelled BRDF data as input for LAI retrieval was explored in this study by coupling the robust PROSAIL model with the hotspot-corrected Ross-Li model using the MODIS BRDF parameter product as input. One of the main contributions of this study is to find a significantly linear relationship between the ALA and f vol of the Ross-Li model in the NIR band, which was used to narrow the search scope of the ALA. The local search method based on the ALA relationship is approximately 11.5 times faster than the wide search method, and therefore, it is promising for rapid large-scale inversion. In addition, the local search method also has a better inversion accuracy of LAI than the wide search method, as validated in Section 4. However, the two search methods also need to be jointly applied, due to the constraint of the f vol -ALA relationship that is applicable for f vol values ranging from 0 to 0.3813. The LAIs subsequently retrieved by combining the two search methods attain good accuracy, especially for high-quality MODIS BRDF data, which presents the importance of quality in modelling multi-angular data. A higher accuracy was achieved for crops than does the MODIS LAI product generated based on limited observations [55], which demonstrates the significance of modelled, sufficient BRDF information in parameter retrieval. For mediumresolution, multi-angular satellites with a few sun-viewing geometries, this novel method can help surface inversion and even attain a higher accuracy. For high-resolution satellites that usually do not have multi-angular observations, these modelled, sufficient BRDF data from the Ross-Li model are promising to act as the compensation.
This study linked the semi-empirical Ross-Li model and the PROSAIL physical model, and the existing model inconsistencies between them may lead to some confusion. The initial two primary components of the Ross-Li model include the volumetric scattering and geometric-optical scattering, which were divided into three parts after employing many simplifications to keep a linear function that is applicable for different spatial scaling [13]. The items associated with background reflectance and a series of parameters of elements (e.g., the height, length and LAI of protrusions) in the whole scene are integrated into three kernel coefficients. Therefore, these non-negative kernel coefficients (i.e., f iso , f vol , and f geo ) are strongly associated with canopy structure, which can roughly indicate the weight of the scattering type. Through numerous validations, the RTLSR model was demonstrated to be robust among a series of kernel combinations, and then it was used as the operational algorithm for generating the global MODIS BRDF parameter product. The hotspot revised version (i.e., the RTLSR_C model) used in this study has significantly improved the accuracy in modelling hotspot reflectance, which also presents a high inheritance from the RTLSR model in the remaining geometries. For the PROSAIL model, more elaborate descriptions of light scattering have been considered, following radiative transfer theory [43,44]. Despite model differences, a previous study represented an overall good BRDF agreement in linking the two models with a wide change in the PROSAIL input parameters [42], which lays a foundation for this study.
In terms of the newly built empirical f vol -ALA relationship, several associated model inconsistencies were discussed here in detail. First, the relationship does not seem to be rigorous, without a clear physical meaning. In fact, both the volumetric scattering Rossthick kernel and the SAIL model are developed based on radiative transfer theory [13,44], which is highly sensitive to the leaf structure inside the canopy, such as the ALA and LAI. Therefore, it is generally understandable for the empirical relationship between the ALA and the f vol that acts as the weight coefficients of the volumetric scattering Rossthick kernel. Second, only single scattering is considered in the Rossthick model, while both single scattering and multiple scattering are included in the SAIL model. Towards the initial volumetric scattering component of the Ross-Li model, the isotropic contribution resulting from multiple scattering has been included in the Lambertian term [13], so multiple scattering has indirectly been considered. Finally, the Rossthick model has also been simplified to describe a spherical distribution of the leaf inclination angle with a large LAI [13], while the PROSAIL model has a wide change in the ALA and LAI. In fact, the broad adaptability of the RTLSR model has been presented in various land cover types from soil to forest-shrub-grass areas [21,57,59,60], and therefore, the model is expected to be applied much more than a specific range of the ALA and LAI, such as the global MODIS BRDF data. This performance may be attributed to the semi-empirical property of the Ross-Li model, where the final fitting accuracies are greatly affected by the observations apart from the theoretical variance of these kernels.
The possible limitations of this study and future expectations of similar studies were also analyzed. First, the following points need to be noted in terms of algorithm accuracy: (1) In the validation of five mainstream global LAI products based on the field-measured crop LAIs used in this study [55], the RMSEs change from 1.24 to 1.77, except for that of the GLASS LAI (i.e., 0.92) for pixel-to-pixel comparisons. In this case, the RMSE of 1.13 for this proposed method shows a relatively good accuracy compared to these LAI products. (2) This study focuses on the improvement of LAI retrieval by using sufficient multiangle reflectances modelled from MODIS BRDF parameters, which was demonstrated by the better estimation accuracy of the novel method than that of the Collection 6 MODIS LAI product generated from observations at a few angles (i.e., RMSE = 1.50). (3) Even so, the accuracy of this study still seems to be large considering the target accuracy (±0.5) required by the Global Climate Observing System (GCOS) [8]. Based on global sites over forest-shrub-grass surfaces, validations of the Collection 6 MODIS LAI product show a smaller RMSE of 0.66 [8]; therefore, the accuracies of this novel method at sufficient sites and vegetation types need to be further investigated. Notably, the reliability of the LAI indirectly validated the retrieval accuracy of the ALA. Then, the potential issue of the Ross-Li model regarding large angles (e.g., SZA values > 70 • ) [71,72] may lead to insufficiency constraints of the corresponding BRDF data when performing retrieval. In addition, the PROSAIL model is more applicable for homogeneous continuous surfaces, such as the crops associated in this study, and geometric-optical physical models (e.g., 5-scale) can be used for heterogeneous discrete canopies, such as forests, in the future. More importantly, the general flowchart of LAI estimation in this study can be applied to different physical and semi-empirical/empirical models [73][74][75][76][77][78][79][80] to retrieve various vegetation parameters, which is one of the main implications for the field of surface parameter estimation. Before model linking, the BRDF consistency should be investigated to reduce model coupling errors [42].

Conclusions
Satellite data play an essential role in surface parameter inversion for monitoring changes in global climate and ecology. Currently, satellite-observed reflectances at a few viewing angles are often directly used for vegetation structure parameter retrieval; however, noise in observations as well as these insufficient sun-viewing geometries always bring difficulties to accurate inversion, generally from physical BRDF models. To address these limitations, modelled sufficient multi-angular reflectances through BRDF models have been widely used as supplemented anisotropy information, which has achieved good performance in the retrieval of several essential vegetation parameters. In this study, we attempted to apply this reconstructed BRDF information to better retrieve LAI from the physical PROSAIL BRDF model using the widely used semi-empirical Ross-Li BRDF model by inputting the MODIS BRDF parameter product. As a contrast, we compared the estimated LAIs with the MODIS LAI product that are retrieved from satellite-observed reflectances at a few angles, which aims to evaluate the potential of the sufficiently modelled BRDF information in the improvement of LAI retrieval. Specifically, an empirical linear relationship between the ALA and f vol of the Ross-Li model in the NIR band was built, which was used for a local search compared to the wide search strategy to determine the optimal combinations of the PROSAIL model parameters. By combining the two search strategies, the novel method attains a good LAI retrieval accuracy, especially given highquality MODIS BRDF data. Both the validation results at the field-measured plot level (i.e., RMSE = 1.28, bias = −0.46, RRMSE = 39.5%) and 30-m LAI maps (i.e., RMSE = 1.13, bias = −0.19, RRMSE = 36.8%) for crops achieve higher accuracies than the MODIS LAI product, which is generated based on observed reflectance data at a few angles.
The good performance of the novel method demonstrates the significance of the modelled, sufficient BRDF information in parameter retrieval, and the obtained ALAf vol relationship can significantly accelerate the inversion process and improve accuracy. Despite several limitations in our method, as introduced in the Discussion section, the developed method is promising for LAI estimation in large homogeneous regions, and the flowchart of this study can be applied to a series of empirical/semi-empirical and physical BRDF models for retrieving various surface parameters. Meanwhile, more algorithm validations at widespread sites for homogeneous canopies need to be performed in future work. In addition, efforts should be made to reduce the uncertainties of the Ross-Li BRDF model at large angles to provide more accurate anisotropy reflectance reconstructions, and geometric-optical physical BRDF models can be joined to retrieve vegetation parameters from heterogeneous surfaces.