Non-Parametric Statistical Approaches for Leaf Area Index Estimation from Sentinel-2 Data: A Multi-Crop Assessment

: The leaf area index (LAI) is a key biophysical variable for agroecosystem monitoring, as well as a relevant state variable in crop modelling. For this reason, temporal and spatial determination of LAI are required to improve the understanding of several land surface processes related to vegetation dynamics and crop growth. Despite the large number of retrieved LAI products and the efforts to develop new and updated algorithms for LAI estimation, the available products are not yet capable of capturing site-speciﬁc variability, as requested in many agricultural applications. The objective of this study was to evaluate the potential of non-parametric approaches for multi-temporal LAI retrieval by Sentinel-2 multispectral data, in comparison with a VI-based parametric approach. For this purpose, we built a large database combining a multispectral satellite data set and ground LAI measurements collected over two growing seasons (2018 and 2019), including three crops (i.e., winter wheat, maize, and alfalfa) characterized by different growing cycles and canopy structures, and considering different agronomic conditions (i.e., at three farms in three different sites). The accuracy of parametric and non-parametric methods for LAI estimation was assessed by cross-validation (CV) at both the pixel and ﬁeld levels over mixed-crop (MC) and crop-speciﬁc (CS) data sets. Overall, the non-parametric approach showed a higher accuracy of prediction at pixel level than parametric methods, and it was also observed that Gaussian Process Regression (GPR) did not provide any signiﬁcant difference ( p -value > 0.05) between the predicted values of LAI in the MC and CS data sets, regardless of the crop. Indeed, GPR at the ﬁeld level showed a cross-validated coefﬁcient of determination (R 2CV ) higher than 0.80 for all three crops.


Introduction
With the global increase in food demand, obtaining timely information regarding crop growth and retrieving detailed data of crop health have become essential aspects for developing strategic food policies and ensuring sustainable agroecosystem management. Reliable monitoring of crop yields at the regional scale can support policy makers in quantifying food supply (GEOGLAM initiative, https://www.earthobservations.org/geoglam.php, accessed on 20 April 2021), while mapping crop conditions at a field scale can assist farmers in agroecosystem management [1]. The leaf area index (LAI), defined as the total one-sided Remote Sens. 2021, 13, 2841 3 of 25 RTM. Moreover, García-Haro et al. [29] used an AVHRR sensor to obtain an operational LAI product at 1.1 km of spatial resolution and 10 days of temporal resolution using GPR. Baret et al. [30] used a SPOT/VEGETATION sensor and, based on a hybrid non-parametric approach (Neural Network and GPR), retrieved a global LAI with 10 days of temporal resolution and 1.5 km of spatial resolution. Despite the large number of retrieved LAI products and the efforts to develop new and updated algorithms for LAI estimation, the available products are not yet capable of capturing the site-specific variability required in many agricultural applications, Therefore, specific LAI data sets with higher spatial and temporal consistency are required. To cope with the spatial and temporal variability of heterogeneous agricultural systems, high spatial resolution satellite systems (10-30 m), such as Sentinel-2, have been gathering growing interest in the field of agronomic research [24,26]. Thus, studies devoted to improving the accuracy and spatial resolution of LAI estimation are needed, in order to respond to the increasing demand for tools to support the site-specific management of crops and landscape [31]. Therefore, this study was focused on the robustness of non-parametric approaches in relation to different sources of variability such as crop species-, growth stage-, farm-, and year-specific conditions, in order to capture site-specific variability.
The objectives of this study were: (i) to evaluate the potential of non-parametric algorithms at pixel level (within field) and at field level for multi-crop and multi-temporal LAI retrieval; and (ii) to test the temporal consistency of the retrieved LAI at field level for crop monitoring over the entire growing season. Thus, to achieve these objectives, we performed an intensive field campaign, contemporary to S2 data acquisition, to collect ground-LAI measurements collected in Tuscany (Central Italy) over two growing seasons (2018 and 2019), including three crops (i.e., winter wheat, maize, and alfalfa), characterized by different growing periods and canopy structures, and considering different agronomic conditions (i.e., three farms in three different sites). Indeed, the final database (ground-LAI and S2 spectral data) is a potential contribution for other studies, thanks to the spatialspectral-temporal characteristics of LAI data.

Materials and Methods
The methodology followed for the assessment of the performances of parametric and non-parametric approaches consisted of four steps: (i) data collection and pre-processing (Section 2.2); (ii) model definition parametric and non-parametric LAI retrieval (Section 2.3); (iii) accuracy assessment (Section 2.4); and (iv) temporal consistency analysis of LAI retrieval (Section 2.5). Concerning MLRAs, further analysis was performed, in order to evaluate the selection of bands in the retrieval process (Section 2.3.2.1). Indeed, the strength of MLRAs, with respect to regression based on Vis, involves the exploitation of the full spectral information of optical data. The spectral band contribution is important to assess such a contribution and to interpret the model results. The flowchart is presented in Figure 1.
operational LAI product at 1.1 km of spatial resolution and 10 days of temporal resolution using GPR. Baret et al. [30] used a SPOT/VEGETATION sensor and, based on a hybrid non-parametric approach (Neural Network and GPR), retrieved a global LAI with 10 days of temporal resolution and 1.5 km of spatial resolution. Despite the large number of retrieved LAI products and the efforts to develop new and updated algorithms for LAI estimation, the available products are not yet capable of capturing the site-specific variability required in many agricultural applications, Therefore, specific LAI data sets with higher spatial and temporal consistency are required. To cope with the spatial and temporal variability of heterogeneous agricultural systems, high spatial resolution satellite systems (10-30 m), such as Sentinel-2, have been gathering growing interest in the field of agronomic research [24,26]. Thus, studies devoted to improving the accuracy and spatial resolution of LAI estimation are needed, in order to respond to the increasing demand for tools to support the site-specific management of crops and landscape [31]. Therefore, this study was focused on the robustness of non-parametric approaches in relation to different sources of variability such as crop species-, growth stage-, farm-, and year-specific conditions, in order to capture site-specific variability.
The objectives of this study were: (i) to evaluate the potential of non-parametric algorithms at pixel level (within field) and at field level for multi-crop and multi-temporal LAI retrieval; and (ii) to test the temporal consistency of the retrieved LAI at field level for crop monitoring over the entire growing season. Thus, to achieve these objectives, we performed an intensive field campaign, contemporary to S2 data acquisition, to collect ground-LAI measurements collected in Tuscany (Central Italy) over two growing seasons (2018 and 2019), including three crops (i.e., winter wheat, maize, and alfalfa), characterized by different growing periods and canopy structures, and considering different agronomic conditions (i.e., three farms in three different sites). Indeed, the final database (ground-LAI and S2 spectral data) is a potential contribution for other studies, thanks to the spatial-spectral-temporal characteristics of LAI data.

Materials and Methods
The methodology followed for the assessment of the performances of parametric and non-parametric approaches consisted of four steps: (i) data collection and pre-processing (Section 2.2); (ii) model definition parametric and non-parametric LAI retrieval (Section. 2.3); (iii) accuracy assessment (Section 2.4); and (iv) temporal consistency analysis of LAI retrieval (Section 2.5). Concerning MLRAs, further analysis was performed, in order to evaluate the selection of bands in the retrieval process (Section 2.3.2.1). Indeed, the strength of MLRAs, with respect to regression based on Vis, involves the exploitation of the full spectral information of optical data. The spectral band contribution is important to assess such a contribution and to interpret the model results. The flowchart is presented in Figure 1. Step 1 (Yellow boxes) data preparation; Step 2 (Blue boxes) model implementation; Step 3 (Green boxes) model accuracy assessment; and Step 4 (Grey boxes) LAI estimation analysis.

Study Area
The study area is located 10 km from Pisa, Tuscany, Central Italy ( Figure 2a). The area is flat and the climate is Mediterranean, with a mean annual precipitation of 907 mm and a mean annual temperature of 15 • C (long-term average 1986-2016). According to land-cover spatial information of the Tuscany regional authorities (http://dati.toscana.it/, accessed on 20 April 2021), in 2017-2019, the three prevalent arable crop categories were: (i) cool-season cereals; (ii) perennial meadows; and (iii) warm-season cereals. In order to cope with the main categories of the area, the following crops were considered in this study: C1-winter wheat (Triticum aestivum L.); C2-maize (Zea mays L.); and C3-alfalfa (Medicago sativa L.). The three crops were monitored along their growing seasons at three different test sites. The three sites were located in three different Farms ( Step 1 (Yellow boxes) data preparation; Step 2 (Blue boxes) model implementation; Step 3 (Green boxes) model accuracy assessment; and Step 4 (Grey boxes) LAI estimation analysis.

Study Area
The study area is located 10 km from Pisa, Tuscany, Central Italy ( Figure 2a). The area is flat and the climate is Mediterranean, with a mean annual precipitation of 907 mm and a mean annual temperature of 15 °C (long-term average 1986-2016). According to land-cover spatial information of the Tuscany regional authorities (http://dati.toscana.it/, accessed on 20 April 2021), in 2017-2019, the three prevalent arable crop categories were: (i) cool-season cereals; (ii) perennial meadows; and (iii) warm-season cereals. In order to cope with the main categories of the area, the following crops were considered in this study: C1-winter wheat (Triticum aestivum L.); C2-maize (Zea mays L.); and C3-alfalfa (Medicago sativa L.). The three crops were monitored along their growing seasons at three different test sites. The three sites were located in three different Farms (

In Situ Measurements
At each site, ground-LAI measurements were collected using the Validation of Land European Remote Sensing Instruments (VALERI) sampling strategy [32,33]. VALERI is based on an Elementary Sampling Units (ESU) upscaling approach, in order to capture the variability between the investigated fields and within the field of each crop, while including the variability within a theoretical pixel [33]. The corresponding scale of ESU is a high spatial resolution pixel of S2, and ESU is defined as a plot of 20 m × 20 m, where an investigated area of 1 m 2 was selected. The ground-LAI measurements were conducted from March 2018 to October 2019. In particular, in the 2018 field campaign, measurements were carried out monthly and in each farm; while, in the 2019 field campaign, weekly

In Situ Measurements
At each site, ground-LAI measurements were collected using the Validation of Land European Remote Sensing Instruments (VALERI) sampling strategy [32,33]. VALERI is based on an Elementary Sampling Units (ESU) upscaling approach, in order to capture the variability between the investigated fields and within the field of each crop, while including the variability within a theoretical pixel [33]. The corresponding scale of ESU is a high spatial resolution pixel of S2, and ESU is defined as a plot of 20 m × 20 m, where an investigated area of 1 m 2 was selected. The ground-LAI measurements were conducted from March 2018 to October 2019. In particular, in the 2018 field campaign, measurements were carried out monthly and in each farm; while, in the 2019 field campaign, weekly measurements were carried out only in F3. On each sampling date, for each field, we collected 12 ESU; where, within each ESU, four (2018) and one (2019) replicate measurements of ground-LAI were sampled. Subsequently, in order to obtain a unique representative value for ESU, the ground-LAI mean was calculated within each ESU. A Remote Sens. 2021, 13, 2841 5 of 25 total of 598 measured ESUs were obtained over the study area. Ground-LAI measurements were performed using SunScan (Delta-T Devices, Cambridge, UK), assuming a random distribution of the foliage (effective ground-LAI) under clear-sky conditions. Moreover, to include bare soil conditions [34], 96 bare soil (LAI = 0.0) ESUs were included in the ground-LAI database (12 for each field, during the germination stage of winter wheat and maize for each sampling year). Therefore, we obtained a database consisting of 694 observations. By means of phenological observations and expert knowledge, we conducted a preliminary screening of ground-LAI anomalies, in order to exclude sampling disturbances arising out of instrument calibration. Moreover, Grubbs' test was performed on ground-LAI measurements data, in order to identify and flag outliers [35]. Both screenings were conducted per field per date, in order to avoid possible sources of disturbance on the day of sampling. Grubbs' test was performed with the "outliers" library in the R environment [36]. After the screenings, a total of 558 representative ESUs were maintained, while 232 ESUs were excluded from the measured ground-LAI database ( Table 1). Beside LAI measurements, the phenological stages of the crops were recorded during each sampling date, based on the Biologische Bundesanstalt, Bundessortenamt, and Chemical Industry scale (BBCH), a German scale used to identify the phenological development stages of a plant [37]. Data of crop phenology were arranged, in order to refer LAI measurements to the main crop development stages. Thus, for winter wheat and maize, four main stages were identified: GE, from germination to full emergence; SE, from early leaf development to complete stem elongation; Fl, from initial stages of flower differentiation to end of anthesis; and FD, from initial stages of fruit development to complete fruit maturity. Conversely, in alfalfa, only two stages were considered: Vg, the vegetative stage, occurring after cut or overwintering; and Fl, from initial stages of flowering to mowing.

Sentinel-2 Data
Copernicus Sentinel-2 (S2) is a satellite mission carrying the MSI multispectral sensor, which is characterised by high spatial resolution (10 m, 20 m, and 60 m), high revisit capability (5 days with two satellites), and a moderately large band set (13 spectral bands) from the visible to short-wave infrared (Table A1) [38,39]. The S2 Level 2A (L2A) images were downloaded from the Theia Land Data Centre, which provides time-series of topcanopy surface reflectance which is orthorectified, terrain-flattened, and atmospherically corrected using the MACCS-ATCOR Joint Algorithm (MAJA) [40]. A total of 37 cloud-free images, collected in correspondence with the in situ monitoring period, were used to analyze the relationship between measured ground-LAI and S2 data. Moreover, the most Remote Sens. 2021, 13, 2841 6 of 25 commonly-used S2 bands for vegetation studies at 10 m (B02, B03, B04, and B08) and 20 m (B05, B06, B07, B08A, B11, and B12) were selected for the analysis (Table A1) [31]. Then, the 10 m spatial resolution bands were resampled to 20 m spatial resolution.
In order to couple the ground-LAI values and the S2 spectral information, the centroid of each ESU was used to extract zonal statistics of S2 time-series using the R software package ''raster" [36,41]. S2 reflectance data values within ±5 days from ground data collection were associated to ground-LAI values. The association of ground-LAI values with the corresponding S2 data was carried out using the SQL database software PostgreSQL 9.5, by joining the two data sets. As a result, a complete SQL database, consisting of 558 records of coupled ground-LAI and S2 data, was obtained for the three crops of each farm in the reference period (October 2017 to October 2019; Figure A1 (Appendix A)).

Parametric Methods
By means of a parametric approach, the empirical relationship between vegetation indices (VIs) and ground-LAI was analyzed, including different sources of variability (crop type, farm conditions, and crop growth). The VIs given in Table 2 were selected, based on previous studies which evaluated visible, red edge, and shortwave infrared as the most effective wavelengths for LAI estimation over different crop types [34,42]. Table 2. Vegetation Indices (VIs).

VIs
Name Formula References Normalized Burn Ratio (B08-B12)/(B08 + B12) [44] Linear and non-linear functions were used for LAI prediction. The linear model (LM) was selected, according to a previous study which demonstrated the good performance of LM for LAI prediction in a mixed-crop scenario [34]. The non-linear model was selected, considering the results of previous studies evaluating the relationship between the ground-LAI and VIs [42,45]. Specifically, the dose-response three-parameter logistic model was selected, in order to evaluate the logistic relationship between ground-LAI and VIs [21]. For the linear model (LM), we adopted the general function of Equation (1) and its inverse (Equation (2)), where a is the intercept and b is the slope: The inverse non-linear function (LogIF d ) of the three-parameter logistic model (Equation (3)) was computed based on the parameterization carried out by the doseresponse "drc" R package [46]. In the equation, d is the curve plateau, which indicates the level at which the VIs saturate; b is the relative slope, which indicates the steepness of the curve; and e is the inflection point, which indicates the point at which the VI value is halfway to its saturation level. Based on Equation (3), the inverse of the three-parameter logistic model was calculated using Equation (4).
Considering that the domain of the inverse function is VI ≤ d, the LAI results are unpredictable above the asymptote d, as the function is not real. In order to overcome the problem of losing information due to the saturation of VIs, we evaluated the option Remote Sens. 2021, 13, 2841 7 of 25 of assigning a default value to pixels with VI > d, which is equal to the maximum value of VI measured below d; while, for VI ≤ d, the pixel value was computed based on the inverse function.
Moreover, to evaluate the robustness of the LM and LogIF d over different sources of variability, the functions were parametrized according to the mixed-crop (MC) data set and the crop-specific (CS) subset. Regression analysis was performed using the "drc" package [47,48] in the R environment [36]. In particular, the general model fitting function dose-response model (drm) was used to fit the regression models and, according to the same parameterization, the inverse functions were calculated.

Non-Parametric Methods
The MLRAs were used for multi-crop and multi-temporal LAI estimation by a nonparametric approach. The architecture of these MLRAs and their training processes were based on iterative regression process, and the core of each algorithm is reported in Table 3. More specifically, GPR provides the predictive mean, as well as predictive variance, maximizing the marginal likelihood in the training set, which is learned by hyperparameters through an appropriate kernel function [49]. BAGTREE builds multiple decision trees by iteratively replacing resampled training data and voting for the decision trees, thus leading to a consensus prediction [50]. BOOST incrementally builds an ensemble by training each new instance to emphasize the training instances which were previously mismodelled [51]. Moreover, to evaluate the MLRA performances over different sources of variability, the algorithms were trained over the mixed-crop (MC) and crop-specific (CS) data set collected during the 2018 and 2019 field campaigns. The MLRAs were run by means of the MLRA toolbox [52] of the Automated Radiative Transfer Models Operator (ARTMO) software [53]. Crop development stages strongly influence the reflectance information provided by RS sensors. Therefore, information on spectral relevance may be a useful tool for understanding several crop-related reflectance conditions. In this study, the spectral relevance of each S2 band (Section 2.2.2) used for LAI prediction was evaluated, according to the crop type and crop development stage. The GPR models were trained over each growth stage within the crop-specific data sets (CS) and over the crops on the mixed-crop data set (MC). We exploited the property of GPR to evaluate the predictive capacity of each single band by a covariance function, defined as follows: Information of spectral relevance of each S2 band was obtained by the hyperparameter sigma (σ), which is the weight assigned to the band (b). The σ b value was provided by the GPR within the MLRA toolbox implemented in the ARTMO software [52]. The toolbox provides an absolute σ b value for each used band. In the supplied hyperparameter, the higher the value of σ b , the lower the relevance of the band. In this study, in order to obtain a positive representation of the spectral relevance, we converted the lower values of σ b into higher ones and calculated the relative σ b . This is detailed in Equation (6), where σ b is the Remote Sens. 2021, 13, 2841 8 of 25 spectral relevance, maxσ b is the maximum value of σ b among the bands used for the GRP training, and sumσ b is the sum of σ b among the bands used for the GRP training:

Accuracy Assessment
The accuracy of parametric and non-parametric methods was assessed at pixel and field level using K-fold cross-validation. The K-fold process consists of dividing the data set into k mutually exclusive groups following a k-fold cross-validation partitioning design [54]. In our case, data were randomly split into k = 3 subsets of equal size, of which, iteratively, two were used for calibration and one for validation. Using this approach, the dependence on a single random partition into calibration and validation data sets was reduced and all observations were used for both training and validation, with each observation used for validation just one time [55]. The validation at field level was conducted using linear regression between ground-LAI values and predicted LAI, averaged per crop, per field, and per growing stage, while the validation at pixel level was conducted at the ESU scale. The coefficient of determination (R 2 ) and root mean square error (RMSE) were calculated to assess the prediction accuracy. The parametric and non-parametric models were selected as those robust and efficient for LAI prediction, considering the average value of RMSE from the cross-validation (RMSE CV ) and the average of coefficient of determination estimated between ground-LAI and predicted LAI (R 2 CV ). Finally, the Wilcoxon signed-rank test was used to compare the distributions of predicted LAI values, according to the crop-specific (CS) and mixed-crop (MC) parameterization, and to identify statistical differences between the two parameterization approaches.

Temporal Consistency
For assessment of the temporal consistency of the estimated LAI against the measured LAI, the 2019 ground-LAI data set was used. First, the LAI time-series from the Sentinel-2 L2A data from October 2018 to October 2019 were retrieved, in order to cover the entire reference period. Then, the temporal consistency was evaluated.
To assess the temporal consistency of the predicted LAI for the temporal profiles, one representative ESU per crop was randomly selected. Subsequently, in order to obtain a full time-series, the ESU profile was interpolated and smoothed using the GPR method. Then, profiles of the retrieved LAI were compared with the measured ground-LAI and the related standard errors. The temporal profile was filtered using the Decomposition and Analysis of Time Series software (DATimeS) [56].

Parametric Model
To evaluate the predictability of LAI by VIs (NDVI, SeLI, and NBR), the accuracy of the three-parameter logistic (LogIF d ) and linear (LM) functions were assessed over the crop-specific (CS) and mixed-crop (MC) parameterizations. Moreover, the LM and LogIF d MC parameterizations were assessed and compared at pixel and field level, in order to assess the suitability of VIs for the prediction of LAI in a mixed-crop scenario. Table 4 shows the results of the accuracy of the parametric functions at pixel level. In general, it was observed that the LM showed statistical differences between the predicted LAI values for both mixed-crop and crop-specific parameterizations (p-value < 0.05). In contrast, LogIF d showed that there were no statistical differences between the two parameterizations based on NBR for winter wheat (C1), NDVI for maize (C2), and SeLI for alfalfa (C3) (p-value > 0.05). Overall, the two parametric methods showed R 2 CV values lower than 0.73 for the three VIs.   Conducting a comparison among all three of the crops, LM exhibited the weakest accuracy in C3, considering all the VIs, with a coefficient of regression lower than 0.17. Analyzing the slope and the intercept of the regression equation, it can be argued that the LM led to overestimation, particularly in the case of low LAI values, independent of the VIs and the parameterization data set used. On the other hand, LogIF d strongly overpredicted when the LAI value was in the high-medium range.
In At the field level, both parameterizations of LM ( Figure A5) and LogIF d (Figure 4) showed an R 2 value higher than 0.8 for both C1 and C2. When LM was used, in both C1 and C2, the highest accuracy for the MC parameterization was obtained with SeLI, whose R 2 CV was equal to 0.88 and 0.82 for C1 and C2, respectively. Otherwise, for C1 and C2 in the MC parameterization, LogIF d showed a higher performance using NBR, presenting R 2 CV values equal to 0.87 and 0.85, respectively.

Non-Parametric Model
Even for non-parametric models, MLRAs were trained over the crop-specific (CS) and mixed-crop (MC) data sets, in order to evaluate and compare the accuracy of LAI prediction at pixel level. Table 5 shows the accuracy metrics of MLRAs at the pixel level. In general, the GPR showed the highest R 2 CV and the lowest RMSE CV for each training data set of all the three crops, except for C2. In fact, when BAGTREE was used in the mixed-crop data set of C2, the R 2 CV was 0.62, compared to 0.59 under GPR, and the RMSE CV was 1.02 vs. 1.06 of GPR. When GPR was used in C3, the estimated values were slightly overestimated; indeed, the values of the intercept higher than 1 confirmed this overprediction. However, despite the overestimation of low-LAI values with GPR, the slope values close to 1 highlighted that there was a positive linear relationship between the measured and predicted values.
Overall, the MLRAs showed a higher accuracy of prediction at pixel level than parametric methods, and it was also observed that GPR did not provide any significant difference between the predicted LAI values according to the MC and CS data sets, regardless of the considered crop. The higher efficacy of MRLAs, in terms of the estimation of LAI under both the MC and CS data sets, was also observed at the field level ( Figure 5). Indeed, GPR obtained an R 2 CV value higher than 0.80 for all the crops.

GPR Spectral Band Relevance
The graphs shown in Figure 6 synthesize the analysis performed to evaluate the relevance of S2 spectral bands to the different GPR models for LAI estimation, according to crop development stages for the "crop-specific model" and crop type for the "mixedcrop model".

Non-Parametric Model
Even for non-parametric models, MLRAs were trained over the crop-specific (CS) and mixed-crop (MC) data sets, in order to evaluate and compare the accuracy of LAI prediction at pixel level. Table 5 shows the accuracy metrics of MLRAs at the pixel level. In general, the GPR showed the highest R 2 CV and the lowest RMSECV for each training data set of all the three crops, except for C2. In fact, when BAGTREE was used in the mixedcrop data set of C2, the R 2 CV was 0.62, compared to 0.59 under GPR, and the RMSECV was 1.02 vs. 1.06 of GPR. When GPR was used in C3, the estimated values were slightly overestimated; indeed, the values of the intercept higher than 1 confirmed this overprediction. However, despite the overestimation of low-LAI values with GPR, the slope values close to 1 highlighted that there was a positive linear relationship between the measured and predicted values.
Overall, the MLRAs showed a higher accuracy of prediction at pixel level than parametric methods, and it was also observed that GPR did not provide any significant difference between the predicted LAI values according to the MC and CS data sets, regardless of the considered crop. The higher efficacy of MRLAs, in terms of the estimation of LAI under both the MC and CS data sets, was also observed at the field level ( Figure 5). Indeed, GPR obtained an R 2 CV value higher than 0.80 for all the crops.  Table 5. Cross-validation results of LAI estimation at pixel level from GPR, BOOST, and BAGTREE. The table reports, for each crop (C1 = winter wheat, C2 = maize, and C3 = alfalfa), the metrics obtained according to the parameterization made on the crop-specific (CS) and the mixed-crop (MC) data sets. The table reports the mean of coefficient of determination (R 2 CV ) and root mean square error (RMSE CV ) estimated from the k-fold cross validation procedure (k = 3); as well as the slope and the intercept of regression lines. p-values < 0.05 (*), < 0.01 (**), and < 0.001 (***) indicate significant differences between the predicted LAI with MC and CS data sets.   Regarding the crop-specific model, when C1 (winter wheat) (Figure 6a) was in the stem elongation (SE) stage, the highest relevant bands were in the red edge (RE) region, as B5 and B6 showed a spectral contribution of 22% and 23%, respectively, for a total of more than 45%. Conversely, during the flowering (Fl) stage, a slight decrease of RE contribution (15% B5 and 15% B6), and a significant increase in the near-infrared (NIR) region (B7, B8, and B8A), higher than 15%, was observed. When C1 reached the fruit development (FD) stage, all bands in the VIS-NIR region had a similar relevance, ranging from 10% to 13% contribution, while those in the SWIR region (B11 and B12) showed values lower than 2%.

Metrics
When C2 (maize) (Figure 6b) was in the SE stage, the highest contribution was observed in the NIR (B7, B8, and B8A) and RE regions (B5 and B6), which reached a relevance higher than 10% (up to 14%), while the lowest contribution was exhibited by the blue band (B2), with 0% relevance. During the Fl stage of C2, B8A (with about 14%) showed the highest percentage of relevance, followed by B7 and B8 (with more than 12%). During the FD stage of C2, B8 showed the highest percentage of relevance (18%), B2 also had a great contribution (17%), as well as B7 and B6 (both about 14%), while the contributions of other bands were all around or less than 5%.

GPR Spectral Band Relevance
The graphs shown in Figure 6 synthesize the analysis performed to evaluate the relevance of S2 spectral bands to the different GPR models for LAI estimation, according to crop development stages for the "crop-specific model" and crop type for the "mixedcrop model". Regarding the crop-specific model, when C1 (winter wheat) (Figure 6a) was in the stem elongation (SE) stage, the highest relevant bands were in the red edge (RE) region, as B5 and B6 showed a spectral contribution of 22% and 23%, respectively, for a total of more than 45%. Conversely, during the flowering (Fl) stage, a slight decrease of RE contribution (15% B5 and 15% B6), and a significant increase in the near-infrared (NIR) region (B7, B8, and B8A), higher than 15%, was observed. When C1 reached the fruit development (FD) stage, all bands in the VIS-NIR region had a similar relevance, ranging from 10% to 13% contribution, while those in the SWIR region (B11 and B12) showed values lower than 2%.
When C2 (maize) (Figure 6b) was in the SE stage, the highest contribution was observed in the NIR (B7, B8, and B8A) and RE regions (B5 and B6), which reached a relevance higher than 10% (up to 14%), while the lowest contribution was exhibited by the blue band (B2), with 0% relevance. During the Fl stage of C2, B8A (with about 14%) showed the highest percentage of relevance, followed by B7 and B8 (with more than 12%). During the FD stage of C2, B8 showed the highest percentage of relevance (18%), B2 also Figure 6. Radar plots of the relative weight (%), estimated from σ values provided by GPR analysis, of Sentinel-2 spectral bands to LAI estimation, according to: (i) different growing stages within the crop-specific data set for C1 (winter wheat), C2 (maize), and alfalfa (C3); and (ii) crop type within the mixed-crop data set. The growing stages for winter wheat and maize are: GE, from germination to emergence; SE, stem elongation; Fl, flowering; and FD, fruit development. For alfalfa, the growing stages were: Vg, vegetative; and Fl, flowering.
Apart from B2 (0%) and B7 (7%), when C3 (alfalfa) (Figure 6c) was in the vegetative stage (Vg), a uniform spectral contribution was observed across all of the spectrum, with a relevance higher than 10%. When C3 was in the Fl stage, the lowest percentage was observed with respect to B4 (6%) and the highest was with B6 (14%).
Analysis of the mixed-crop data set (Figure 6d) by GPR showed that the RE and NIR regions provided the largest contribution (about 13% for the spectral bands B5, B6, B7, B8, and B8A) for C1 LAI estimation. Regarding C2 LAI estimation, the most relevant bands were B6, B7, and B8A, with percentage higher than 13%. In C3, the LAI model exploited all of the spectrum, with a similar contribution (higher than 10%); except for B2, which was not relevant (0%).

Temporal Consistency of LAI estimation
Using the data of the 2019 sampling campaign, the parametric function and nonparametric algorithm with the highest accuracy for LAI retrieval were selected for the temporal assessment. Figure 7 show the dynamics of the predicted LAI values compared with the measured ground-LAI of one representative ESU for each crop by LM, LogIF d , and GPR. In general, the predicted LAI profiles obtained from the GPR exhibited the closest predicted LAI values to the measured ground-LAI, while the parametric approaches based on LM and LogIF d showed different consistency of temporal dynamics, according to the crop type.
x FOR PEER REVIEW 16 of 25

Discussion
Considering spectral bands (e.g., near-infrared), indicators (e.g., NDVI), and final useful products (e.g., land-cover), there is a need to obtain valuable information for the Specifically, it was observed that LogIF d and GPR in the stem elongation (SE) stage of winter wheat (C1) showed a similar behavior, yielding predictions very close to the ground-LAI values. Conversely, the LM exhibited the most distant predictions to the measured LAI. In C1 flowering (Fl) and fruit development (FD) stages, the closest predicted LAI values to ground-LAI were obtained using the GPR.
In maize (C2), GPR showed a LAI dynamic very similar to the one measured in situ, while LM and LogIF d overestimated the LAI, especially during SE. During the C2 FD stage, GPR and LogIF d showed a similar behavior, both underestimating the predicted LAI values, while LM exhibited a slight overestimation tendency. In alfalfa (C3), during the earliest vegetative (Vg) stage, LM showed the closest predicted LAI values to ground-LAI. During the Fl stage, GPR and LM showed similar performance.

Discussion
Considering spectral bands (e.g., near-infrared), indicators (e.g., NDVI), and final useful products (e.g., land-cover), there is a need to obtain valuable information for the development of LAI, as elaborated for the Copernicus indictors and products. Thus, assessment of the statistical relationship between ground-LAI and satellite remote sensing data for LAI prediction requires assessing the accuracy of LAI retrieval methods, according to ground-LAI estimation errors [42,57]. In this study, we measured ground-LAI for three crops, and showed a wide range of values, compared with those reported in the literature. Specifically, Revill et al. [58], using SunScan, exhibited a lower range of ground-LAI values (min, 0.5; max, 3.5) for winter wheat. However, our field data covered the entire cycle, from emergence to maturity, and agreed with Upreti et al. [26] (min, 1; max, 6.5).
Regarding maize, Facchi et al. [59], using LAI-2000, a ceptometer, and a Hemispherical camera in an experiment to compare (with destructive sampling measurements) the range of maize ground-LAI, the results were comparable to those we sampled in the field (min, 1; max, 5). Finally, for alfalfa, Verger et al. [60] measured a range of ground-LAI (min, 0.8; max, 6.5) that was lower with respect to our collected maxima (LAI > 8) using the LAI-2000 instrument. From this literature comparison, the only anomalous apparent difference of ground-LAI values was observed for alfalfa. However, the average value was coherent for different ESU in the same field for the 2018 sampling at flowering date (just before mowing) and, so, it was not in disagreement with the analysis. This difference, with respect to the literature, could arise from the uncertainty due to the optical instrument, the phenological stage, and the crop reflectance response [61,62].
Uncertainty of LAI estimations could arise from many factors, including crop type, farm management, and temporal variability of ground-LAI measurements [42,63].
In this work, we observed that, when using parametric approaches, SeLI proved to be the most suitable VI for LAI prediction under a mixed crop scenario. This result is in close agreement with previous findings highlighting that the combination of NIR and RE may provide more accurate LAI estimates for different crop types [18,34]. However, in the assessment of parametric methods, we also demonstrated that, at pixel level, the accuracy of LAI estimation was strongly affected not only by the VI selection but also (and not in a negligible manner) by the regression function as well as the parameterization data set [20,63]. The results of VIs also demonstrated the better LAI prediction ability of SeLI for wheat and maize, compared to alfalfa. The different LAI prediction ability of Vis, according to crop types, is in agreement with Herrmann et al. [64]. In this regard, Dong et al. [17], using RapidEye reflectance data and ground-LAI data of different crop types, demonstrated that VIs based on visible and RE regions are generally affected by chlorophyll, water content, and the structural properties of leaves and, therefore, the LAI predictability of a VI may vary markedly among crops and growing stages, according to canopy characteristics. Furthermore, when the crop-specific and mixed-crop parameterization were compared under both regression methods based on linear and logistic models (LM and LogIF d ), we demonstrated that the parameterization data set strongly influenced the LAI prediction and its accuracy. This result built on the study of Nguy-Robertson et al. [20] who, using hyperspectral data, found that RE-based VIs are little affected by crop type and, thus, may facilitate the prediction of LAI for different crops characterized by different canopy structure. This finding is very important for the future availability of operational hyperspectral missions, such as the foreseen Copernicus CHIME and NASA SBG [31,65]. Thus, to obtain an accurate LAI prediction, the selection of suitable VIs is critical, as ideal VIs should be sensitive to the ground-LAI but insensitive to interference factors (e.g., the soil background, canopy structure, and chlorophyll content) [66]. Thus, despite the fact that SeLI yielded the most accurate results, it was also demonstrated that, due to the above-mentioned uncertainties, it might be not as accurate at pixel level.
Moreover, we showed how the three MLRAs (and, in particular, GPR) were able to exploit all information available from a multispectral data set, thus providing more accuracy in LAI prediction compared to VIs. Indeed, considering the relative weight of bands in the GPR algorithm, it was observed that, although the NIR and RE bands were the most relevant for LAI estimation, the other bands also contributed-varying according to the crop type and development stage-to LAI prediction. This finding is in agreement with previous outcomes. Delloy et al. [67] showed how the RE S2 bands (B5, B6, and B7) can improve winter wheat LAI estimation by using a hybrid approach (RTM + neural network). Verrelst et al. [13] compared different MLRAs, using simulated Sentinel-2 reflectance data over different crops types, and concluded that GPR was the most effective algorithm for LAI retrieval. However, despite the promising results of MLRAs, it was also observed that the performance was influenced by the training data set (i.e., MC vs. CS). This is in agreement with Mao et al. [25], who tested the influence of sample size on MLRA performances and highlighted how its influences the accuracy of the algorithms.
Regarding the temporal consistency of LAI estimation, our results showed that LAI retrieved by the parametric approach, at pixel level, was less suitable, with respect to the non-parametric approach. Indeed, the parametric method based on VIs showed a low accuracy, in terms of representing variability within the field, due to a saturation effect occurring especially at high vegetation density. In contrast, GPR allowed us to point out such variation, regardless of the crop development stage [58]. The low ability of parametric methods to account for the within variability of LAI was evidenced by the weak metrics of cross-validation carried out at pixel level. In particular, LM made it evident that parametric approaches may lead to the overestimation of LAI at early stages of crop development and underestimation at full canopy development. In this latter condition, VIs showed their limit in detecting LAI variation, due to the well-known issue of saturation [45]-a limit of VIs that was particularly exacerbated in the case of alfalfa, which reaches canopy closure, after resprouting, in less time compared to a winter cereal (wheat) and a row crop (maize). Conversely, GPR showed a high ability to detect LAI variation at the pixel level, regardless of the development stage, vegetation density, and crop type.
It is well-known that the canopy reflectance is affected by several biophysical and biochemical variables and, thus, the regions of the reflectance spectrum can be associated with different vegetation properties [68,69]. Therefore, exploitation of the full spectrum with non-parametric methods can improve the quality of LAI retrieval [13]. Indeed, our results demonstrated that the GPR outperformed the parametric methods; in addition, it was the most accurate MLRA for LAI prediction at both field and pixel level. The results of this study showed that the VI-based parametric method had a lower accuracy for LAI retrieval than MLRAs. These results suggest that GPR based on Sentinel-2 multispectral images is promising for crop monitoring, from a multi-crop mosaic scenario perspective. Further work will involve applying the MLRAs trained in this work to verify the model stability when applied to an independent data set; this analysis will allow for a full assessment of the robustness and exportability of the model developed.
The capacity of MLRAs to deal with full spectral information is a promising aspect that makes these approaches candidates for the investigation of new-generation hyperspectral data available from ASI-PRISMA, as well as those expected from the foreseen Copernicus CHIME and NASA SBG missions.

Conclusions
In the present study, the capability of different parametric and non-parametric methods for the retrieval of LAI for different arable crops using Sentinel-2 data was assessed. The accuracy and robustness of LAI estimates were compared, based on repeated in situ ground-LAI measurements throughout the crop growing season. Regarding the VI-based parametric regressions, the normalized index SeLI was evaluated as more suitable (i.e., than NDVI and NBR) for LAI retrieval at field level, providing good evaluation metrics by the cross-validation analysis for winter wheat and maize. However, VI-based parametric methods were shown: (i) to be unsuitable for LAI retrieval of alfalfa and mixed crop scenario; (ii) to have a very low accuracy for LAI retrieval at pixel level; and (iii) to have an accuracy of prediction the largely depends on VI selection, the fitting function, and the parameterization data set.
Among the non-parametric regression methods evaluated, the best-performing MLRA belonged to the kernel machine learning regression algorithms. Indeed, GPR was evaluated as the best-performing algorithm for LAI prediction, for the three arable crops evaluated. Using GPR, Sentinel-2 imagery can be used to map the spatial variability of the LAI of different arable crops, having prediction accuracy which is very high at the pixel level regardless of the crop type, growth stage, and the training data set. Moreover, GPR analysis of spectral bands in different phenological stages provided information on the relevance of relative bands in contributing to LAI prediction. However, further studies are required to fully assess the potential of GPR across different crops in different areas, as well as under contrasting agronomic conditions.        near regression between ground-LAI and VIs (NBR, NDVI, and SeLI), according to the Linear function (LM). orresponds to the regression curve, while colored dots are the three crops: yellow for winter wheat (C1), blue ), and orange for alfalfa (C3). The shapes represent the farms (F1, F2, and F3).