In this study, two types of noise were tested: multiplicative wavelength-independent (MI) and additive wavelength-independent (AI) noise [61]. MI is dependent on the BRF. Its term is larger for NIR compared to red spectral bands for typical vegetation spectral responses. MI and AI were added to the RTM spectral bands:

$$
\rho' = \rho\_{\rm RTM} + \rho\_{\rm RTM} \cdot \varepsilon\_{MI} + \varepsilon\_{AI} \tag{2}
$$

where *ρ* is the noise contaminated spectral band, *ρRTM* is the RTM spectral band output, *MI* the MI noise term with *MI* ∼ N (0, *σMI*) and *AI* the AI noise term with *AI* ∼ N (0, *σAI*). Apart from noise free, realisations with MI noise of 0.05, 0.1, 0.2 and 0.3, and AI noise of 0.05, 0.1 and 0.2 were tested. The 0.05 noise level was motivated by the Sentinel-2 mission requirement of 5% error on SR [34]. The other noise levels were pessimistic variations. However, as mentioned before, the noise term has multiple purposes, so that the mission requirements can only be an indication.

#### *3.5. Solar Zenith Angle*

Illumination conditions greatly affect the reflectance of canopies [5,62,63]. As the SZA changes over the course of the year, the internal canopy shadowing varies. Furthermore, SZA is an easy to obtain feature, as it is solely a function of location and time. Therefore, SZA was incorporated in the training workflow to test if it improves LAI prediction. SZA was calculated for local overpass times of the respective missions with the R package RAtmosphere ([64], https://cran.r-project.org/web/ packages/RAtmosphere). For the respective realisations, it was treated as an extra training feature next to spectral bands.

#### *3.6. Machine Learning Regression Algorithms*

Studies on MLRA typically test a range of algorithms to explore their respective (dis-)advantages and cross-comparison results. This was adopted in this study as well. All models were trained to predict LAI, while the independent variables depended on the learning realisation. Multi-variate OLS regression was chosen as a benchmark method. For neural networks, the classic MLP was used (e.g., [19]). In particular, this was the implementation of the Stuttgart Neural Network Simulator in the R package RSNNS (https://cran.r-project.org/web/packages/RSNNS). Networks with *n* + 1 neurons in a single hidden layer were trained, where *n* corresponded to the number of independent variables [19]. Random Forest was selected as an RT algorithm [65,66] and used as implemented in the R ranger package (https://cran.r-project.org/web/packages/ranger). The forests were grown with 500 trees.

Furthermore, three kernel-based methods were used. This type of regression methods translates the—possibly nonlinear—regression problem from the parameter space into a higher dimensional feature space, where it can be solved linearly. Kernel functions implement a notion of similarity function. SVR [67], KRR [14] and GPR [68] with Radial Basis Function (RBF) kernels were tested here. The kernel *σ* hyperparameter was estimated with the sigest function from the kernlab package (https://cran.r-project.org/web/packages/kernlab/). Although developed to estimate *σ* for SVR, results in initial tests were promising for KRR and GPR. For further reading on MLRAs in biophysical parameter estimation, the reader is referred to Verrelst et al. [23].

The general workflow for MLRA application typically involves splitting of the feature database into training and validation sets to tune model parameters. Here, five-fold cross-validation was performed during the tuning process. This was based on the training dataset, which held 2500 samples of the RTM-based database. Next, the model performance for the best tuned parameters was evaluated with the test dataset, which held 500 samples of the RTM based database. This set was never seen by the models during training. Finally, the models were applied on the actual sensor observations and compared with the validation dataset (Section 3.7). However, as we inverted observations of three different sensors, in fact, for each realisation, three separate models were trained.

#### *3.7. Ensemble Analysis and Validation*

In order to analyse how well the MLRAs were able to learn the RTM-produced band-LAI relationships, the test error was evaluated. However, this error represents only the theoretical performance in case the RTM produces true results for the scene and the induced noise properties correspond to the noise of the actual spectral observations.

For validation purposes, the advantages of the TLS time series—i.e., high precision due to independence of illumination conditions [69,70]—was combined with the direct estimation of LAI with the litter-traps. Litter-traps are considered among the most accurate methods in terms of absolute LAI for forest canopies [71]. This approach follows the suggestion of Woodgate et al. [72] to calibrate TLS with other techniques. Specifically, the TLS time series was scaled with the litter-trap total LAI separately for each plot:

$$LAI\_i = \frac{PAI\_{TLS}i - \min(PAI\_{TLS})}{\max(PAI\_{TLS})} \cdot LAI\_{LT} \tag{3}$$

where *LAIi* is the LAI at time *i*, *PAITLS* is the TLS derived PAI time series and *LAILT* is the litter-trap LAI. The observations were averaged per plot and linearly interpolated to obtain a continuous time series.

For each realisation, the time series of predicted LAI was compared with this validation time series. The performance metric was the RMSE:

$$RMSE = \sqrt{\frac{1}{t} \sum\_{i=1}^{t} (LAI\_{realisation,i} - LAI\_{valid,i})^2} \tag{4}$$

where *t* is the length of the time series, and *LAIrealisation*,*i* and *LAIvalid*,*i* the realisation and validation LAI at time *i*, respectively.

#### **4. Results and Discussion**

#### *4.1. Validation Time Series*

Figure 6 shows the derived validation time series. The seasonal pattern with a fast spring leaf flush in May and autumn leaf falling in November dominated the temporal behaviour. Calders et al. [35] observed this speed in spring leaf flush for a mixed oak forest in the Netherlands. Maximum LAI of 6.1 m2 m−2 was reached on 26 May in plot A. Plots B and E showed LAI values just below 6.0 m2 m−2. Plot C had a larger gap in its centre so that overall LAI was lower there. Measured LAI for plot D was about two units lower than for the other plots with a peak of 3.2 m2 m−2. This was due to the age composition of this plot, which was dominated by younger trees. The maximum LAI compares well with the results of Leuschner et al. [73], who measured LAI by litter-traps in 23 mature Beech stands in Germany. They found an average LAI of 7.4 m2 m−2 with a range between 5.6 to 9.5 m2 m−2.

Another feature is the slow decrease in LAI starting in August that could be observed in all plots. After the 2016 growing season, few brown leaves were still remaining on the trees until new leaves flushed in 2017. Overall, the obtained time series show the expected dynamic behaviour of the canopy during spring and autumn (Figure 6).

With respect to the uncertainty of the validation data, the standard deviation of the mean for the litter-trap samples was calculated as 0.43, 0.25, 0.41, 0.32 and 0.24 m2 m−2 for plots A to E, respectively. GCOS specified an accuracy of 0.5 m2 m−2 as a target for LAI products for local and regional applications [38]. However, this is the requirement for the final LAI products, so that the achieved uncertainties are rather high for a validation measurement. This affected especially the realisation evaluations in terms of RMSE.

**Figure 6.** Speulderbos LAI time series for 2016 derived from TLS, points are observations. Lines are linear interpolations. Interpolations for outside of the measurement campaign were performed with the values from the previous and next year, which are added to the graph for clarity.

#### *4.2. Leaf Biochemical Parameters' Retrieval*

The leaf chemistry assessment indicated that some leaf components remained stable over the course of the season, while others showed a dynamic behaviour, resulting in multi-modal distributions (Figure 7). The static ones were the leaf structure parameter (N) parameter, the dry matter content *Cm*, and to some extent the equivalent water thickness *Cw*. The mean of the carotenoids remained stable, but its variance was increasing at the last sampling day. *Cab* showed clear dynamics, with a strong decrease during the last sampling day. When compared to the readings of the SPAD meter, *Cab* retrievals showed a quadratic relationship (Figure A1), which is confirmed by other studies [74,75]. This strong relationship supports the validity of the *Cab* retrievals, and thereby the retrieval of the other biochemical constituents since they were inverted simultaneously. On the basis of these results, it was decided to constrain the training with fixed, central values of the N parameter, leaf carotenoid content (Car), *Cw* and *Cm*, but vary *Cab* and leaf brown pigment fraction (*Cbrown*) as given in Table 2.

**Figure 7.** Leaf chemical properties based on PROSPECT inversions of 173 leaf samples. Leaf properties are leaf structure parameter (N), Chlorophyll a and b (*Cab*), leaf carotenoid content (Car), leaf dry matter content (*Cm*), leaf water content (*Cw*), and leaf brown pigment fraction (*Cbrown*). Stacks are coloured by acquisition date. Solid lines are fitted density models. Ordinate axis represents observation counts and modelled counts for the observations and the fitted models, respectively.

#### *4.3. RTM Sensitivity*

PROSAIL's and DART's sensitivity to their input parameters is depicted in Figure 8. In case of PROSAIL, BRFs in the visible bands were primarily driven by *Cab* with a contribution of 74.8 and 49.4% in B03 and B04, respectively. This extended into the first red edge band B05, but strongly decreased to 6.1% at B06. In the red edge and NIR bands, LAI was the most important parameter with a relative contribution of 68.5% for bands B06 to B8A on average. This sensitivity is the reason why LAI retrieval relies on the NIR domain. The SWIR bands were mostly dependent on leaf water content *Cw*, which had 52.8 and 46.2% contribution in B11 and B12, respectively. These results are in line with those of Jacquemoud et al. [5].

On the other hand, DART's output sensitivity was dominated by canopy structural parameters (LAI and crown diameter) in the NIR spectral outputs (Figure 8). The contribution of LAI was maximal in B04 with 64.1%, while that of crown diameter was maximal in B06 with 52.2%. Their combined contribution was minimal in B12 with 12.7%. In contrast to PROSAIL, DART showed some sensitivity towards SZA, which was 5.1% on average in bands B05 to B8A. This reflects the effect of shadowing and DART's vertical heterogeneous character.

**Figure 8.** Global sensitivity of S2A MSI spectral bands to PROSAIL and DART input parameters according to Sobol' Indices. Total Sobol' Indices were normalised per band to sum up to 1. For band specifications, see Table 1 and for parameters see Table 2.

#### *4.4. Impact of Training Scheme Features on Prediction Performance*

This section presents the single inversion workflow features and their role for predicting LAI, starting with a comparison across domains in the following paragraphs. Since the RMSE results were not consistently normal distributed within ensembles, the median was calculated for all realisations that implemented a specific feature. In this way, the feature's role could be evaluated. It should also be noted that RMSE refers to the validation error. Only in some cases was the training error evaluated, but this is always explicitly mentioned.

Table 3 summarises the effects of all features on the validation RMSE. The 25 and 75% Interquartile Range (IQR) is listed as a quantification for the variability within the ensembles. The feature with the strongest influence on performance in terms of RMSE was the adding of AI noise to the RTM generated spectral outputs. Adding 5% AI noise decreased the median RMSE from 2.33 m2 m−2 for no noise to 1.25 m2 m−2. Realisations with 5% AI noise also achieved the overall lowest RMSE. The second-most important feature was the choice of MLRA. Here, realisations varied between RMSE of 2.04 m2 m−2 for MLP, and RMSE of 1.52 m2 m−2 for SVR. Restraining the training database with prior information on leaf biochemical contents generally decreased prediction performance.

The following sections present the performance results of all training features in more detail and elaborate on their first order interactions, whereas typically only the strongest interaction is discussed. Interactions were investigated by comparing the respective groups of realisations that implement the features. For example, if feature A has two realisations A1 and A2, and its strongest first order interaction feature B has B1 and B2, all realisations that implement them (A1/B1, A1/B2, A2/B1, A2/B2) were compared to each other. The strongest interaction (feature B in the previous example) was identified as the one that varies performance in terms of validation RMSE the most after the feature under consideration (feature A in the previous example). For example, when the leaf chemical prior was inspected, the feature that produced the largest differences among the realisations was identified as the strongest interaction, which was the RTM feature.

**Table 3.** Validation median RMSE and Interquartile Range (IQR) (25 to 75%) for training features. ΔRMSE refers to the difference to the best feature per group.

#### 4.4.1. Leaf Biochemical Prior

Among other information, Figure 9 summarises how the leaf biochemical prior information affected the training results. PROSAIL performance was generally more affected by introducing prior information than DART performance (Figure 9): median RMSE decreased from 1.91 to 1.38 m2 m−2 when prior information was included in the test data sets, while it increased from 1.43 to 2.37 m2 m−2 in case of the validation data. The former can be explained by the sensitivity of PROSAIL to biochemical parameters *Cab*, *Cm* and *Cw* (Section 4.3). Using prior information effectively decreases the parameter input space that the MLRA has to learn, thereby making it easier for the MLRA. This also made it possible to reach testing RMSE as low as 0.01 m2 m−2. However, inversion of actual observations was impaired by the constraint of the leaf chemical parameter space. This may be due to the fact that leave chemical properties were not representative for the whole study area, as leaves were only sampled from two beech trees, while oak was also present and environmental conditions may change the chemical composition throughout the study area. Additionally, the way the constraint was implemented—as mean estimates allowing no deviation—may also play a role.

Contrary to this, the DART inversion performance was less sensitive to leaf biochemical parameters (Section 4.3). Median performance in terms of RMSE was similar at 2.29 and 2.16 m2 m−2 for test data and improved slightly with a decrease in RMSE from 1.55 to 1.46 m2 m−2 for the validation data set when introducing leaf chemical information. Thus, reducing the input parameter space had a small positive effect during training and validation.

**Figure 9.** Prediction performance in terms of testing and validation RMSE with dependence on the chosen RTM, AI and MI noise level, and bio-chemical prior. *RMSEmedian*, *IQRRMSE* and *RMSEmin* refer to the validation error in the respective panel cell. Grey line is 1:1 line. In case of the validation results, 14 realisations were trimmed with RMSE larger than 10 m2 m−2 (all with 0 % AI noise) because they prevented proper display.

#### 4.4.2. RTM Choice

Among other features, Figure 9 compares the inverse model realisation in terms of their used RTM. Both testing error, i.e., the error based on RTM produced samples, and validation error, i.e., the error based on the validation time series for the Speulderbos site, are shown. Overall, PROSAIL and DART achieved a median test RMSE of 1.63 and 2.21 m2 m−2, while the validation RMSE was 1.93 and 1.51 m2 m−2, respectively. The lower performance of DART on the testing samples can be explained with its additional freely varying parameter, the crown diameter. However, with this parameter came the capability to model crown gaps (Figure 5), which led to the better performance in terms of validation RMSE.

Apart from the overall better median RMSE of DART, using this RTM generally decreased spread of error for realisations that implemented this RTM, at least when some AI noise larger than 0% was chosen. For example, the validation RMSE IQR for DART with 0% AI noise was 0.23 m2 m−2, while it was 1.20 m2 m−2 for PROSAIL. This means that choosing this DART reduced the importance for a particular choice of the other training features.

There were realisations for which the testing RMSE was disproportionally larger than the validation RMSE. In fact, this was the case for 43.3 and 77.7% of the PROSAIL and DART cases, respectively (Figure 9). Under circumstances where data of the same origin would have been used, this would be unlikely to occur, especially in scenarios with added noise. However, training in this study was based on RTM output and validation on field acquired data.

#### 4.4.3. Noise Scenarios

Both the AI and MI noise had different effects on the testing and validation error. Generally testing errors increased with increasing noise—in the case of AI from 1.29 m2 m−2 for 0 % to 2.46 m2 m−2 median RMSE for 20% AI noise, and in the case of MI from 1.89 m2 m−2 for 0% to 2.09 m2 m−2 median RMSE for 20 % MI noise. This showed the general effect of the noise to blur the relationship between spectral output and associated LAI, and prevent the MLRA to learn the true RTM produced pattern.

AI noise was generally more successful at reducing validation RMSE (Table 3): it decreased median RMSE by 1.08 m2 m−2 (from 2.33 m2 m−2 for 0% AI noise to 1.25 m2 m−2 for 5 % AI noise). MI noise reduced median RMSE only by 0.10 m2 m−2 (from 1.73 m2 m−2 for 0 % MI noise to 1.63 m2 m−2 for 20% MI noise).

Considering AI noise alone, the addition of any in comparison to 0% AI noise strongly changed the distribution of realisations in terms of RMSE (Figure 9). This did not only include the best, but also low performing results. More precisely, AI noise prevented occurrence of very bad results. For the 0% AI noise level, worst performance reached up to 29.29 m2 m−2 and the 95% quantile lay at 10.4 m2 m−2, while for 5% these statistics were 2.66 and 2.48 m2 m−2, respectively. Here, the added noise prevented the MLRA from over-fitting on clean RTM simulations.

Noise was also found an important training workflow element in PROSAIL-Look Up Table (LUT)-based inversions of observations from agricultural crops by Rivera et al. [29] and Verrelst et al. [30]. They added up to 50 and 30% noise to PROSAIL spectra, respectively, but did not specify how the error was implemented. However, the required magnitude of noise corresponds to the MI noise in this study, which was optimal at 20 to 30%. Baret et al. [19] added 0.04 absolute Gaussian white (AI) noise in a global retrieval workflow based on PROSAIL. This is in line with the optimal 5% AI noise in this study. Koetz et al. [28] adopted a wavelength-dependent, relative noise term of maximal 10% in the 444 nm band. Both Koetz et al. [28] and Baret et al. [19] based their choice for a specific noise level on experience of observation errors, but do not evaluate other possibilities.

Baret et al. [19] argued that the quantification of the error term is difficult, as it includes errors stemming from the radiometric calibration of the sensor, Bidirectional Reflectance Distribution Function (BRDF) normalisation, atmospheric correction, cloud residuals and the RTM representativeness for the actual canopy. Additionally, the interaction of these single terms plus the properties of the used MLRA in an inversion workflow complicates the choice based on experience of errors of the sub-systems. Surely, it is more practical to conduct a sensitivity analysis over validation samples rather than characterising the sensor-inversion system in detail. Moreover, noise terms need to be defined properly to compare them across studies.

#### 4.4.4. SZA

As shown in Table 3, realisations that made use of SZA as an extra training feature achieved overall 0.06 m2 m−2 smaller median RMSE than realisations that did not include SZA. A Wilcoxon signed-rank test confirmed that the two groups differed significantly (*p* < 0.01). This difference was most prominent with the multiple MLRAs. For example, SVR benefited strongest with a decrease in median RMSE by 0.07 m2 m−2. This can be explained by the richer feature space that the MLRAs had available for learning. Additionally, SZA correlated with the general phenological patterns of the study area with low SZA in summer.

Strategies to include SZA into inversion workflows of multi-temporal observations are not consolidated yet. Koetz et al. [28] computed independent LUTs for different observation dates and consequently SZAs, but they did not investigate the error that would occur if they would not have done so. However, training separate models for multi-temporal time series with multiple observations per year is undesirable due to computational load and the checks that would be necessary to ensure consistent model properties. Campos-Taberner et al. [16] conducted PROSAIL inversions over a full season of L8 OLI, L7 ETM+ and Satellite Pour l'Observation de la Terre (SPOT) High Resolution Geometric (HRG) sensor data and they mention solar-sensor geometry as inputs for PROSAIL, but do not elaborate how these parameters were dealt with in the training database. Baret et al. [19] included SZA as a training feature in their neural network based inversion for the VEGETATION based CYCLOPES LAI product, but again did not evaluate this strategy. However, as the CYCLOPES product has global extents, it spans several degrees of latitude, leading to different regimes of illumination dynamics over the year. Thus, SZA was given importance in past studies and showed some importance here as well, but has not been yet evaluated for global LAI products. At least adding SZA as an additional training feature is an easy implementable option. In addition, SZA is efficient to compute, as only the location and observation time are needed to calculate it. However, interactions with other parameters that change over the time of the year such as soil background and LAI itself also need to be considered.

#### 4.4.5. MLRA

As demonstrated in the discussion of the other training features, the choice for a specific MLRA was not the most important factor affecting validation performance in this study (Table 3). However, studies employing MLRA typically compare various algorithms (e.g., [14,16]). Figure 10 gives an overview over all realisations grouped by their used MLRA. The MLRAs reached best (and median) RMSE of 0.91 (1.57), 0.92 (1.63), 0.93 (1.52), 0.95 (1.60), 1.02 (1.72) and 1.08 (2.04) in case of RT, KRR, SVR, GPR, OLS and MLP, respectively. Hence, even though RT and KRR produced the best, SVR produced the overall best realisations. Maximum differences among the MLRAs in median RMSE of 0.52 m2 m−2 were found between SVR and MLP realisations. There were also 18 realisations using OLS and seven using KRR that exceeded RMSE of 5 m2 m−2, all of which were training without noise applied to the training spectral features.

Apart from the MLRA validation, performance processing time is an important property especially for routine and large scale production. The time required for the training of the described realisations in this study was on average 0.03, 7.18, 2.29, 3.17, 125.37 and 123.81 s for the OLS, RT, MLP, SVR, KRR and GPR, respectively. This is in contrast to Verrelst et al. [14] who found KRR and GPR required around the tenth of the time of a neural network or an SVR. However, their models were implemented with Matlab, while this study used R. Additionally, the implementations of RT, SVR, KRR and GPR in this study could make use of parallelisation. Concerning time required for prediction of 10,000 random samples, the models needed 0.08, 0.09, 0.09, 0.15, 0.30 and 0.19 s in the case of OLS, RT, MLP, SVR, KRR and GPR, respectively. Hence, implementation details and optimisation can play a significant role in processing time.

**Figure 10.** Violin plots [76] of prediction performance for the different MLRAs. It should be noted that 25 realisations were trimmed with RMSE larger than 5.0 m2 m−2 (18 for OLS and 7 for KRR) because they prevented proper display.

#### *4.5. Best Performing Feature Combination*

Figure 11 shows the best performing realisation that reached RMSE of 0.91 m2 m−2. It was built with an RT based on a DART produced database restricted with prior information from the leaf sampling, 5% AI and 20% MI noise, and with SZA as an additional training feature. Maximum LAI of 5.80 m2 m−2 was reached in plot A on 17 July. LAI before May 1st was on average (0.27 ± 0.31) m2 m−2, and thereafter (4.81 ± 0.59) m2 m−2 until end of October. The predicted LAI explained 84.9% of variation of the reference time series.

In general, these results for a single time series are in the range of previously published results. Schlerf and Atzberger [77] achieved 0.66 m2 m−2 with a two-band combination chosen from simulated Landsat TM bands over a Norway Spruce site. For beech canopies within the same site, Schlerf and Atzberger [78] report 2.12 m2 m−2 RMSE with a multi-spectral, near-nadir viewing set-up. Brown et al. [18] obtained a RMSE of 0.47 m2 m−2 for a beech site in Southern England. All three studies used the INFORM RTM to create the training database [79].

However, as can be observed in Figure 11, the predictions showed different biases for the single plots. In fact, summer LAI was underestimated on average by 0.74 m2 m−2 in plots A, B, C and E, and overestimated by 1.53 m2 m−2 in plot D. The bias in plot D could result from the structure of the plot, which consisted of more young trees compared to the other plots. The position of the litter traps was chosen close to the TLS scan positions, which were possibly not representative for the whole of the plot. Additionally, RMSE as the choice of error metric during retrieval evaluation does not allow the assessment of bias. In fact, additional error metrics would be needed to characterise the bias as well as temporal consistency of the time series.

**Figure 11.** Best performing realisation in terms of RMSE over the 2016 study period. Black solid line is the validation time series.

#### **5. Conclusions**

The Sentinel-2 mission provides a new science-grade data stream for monitoring of dynamic vegetation behaviour. Together with other missions such as Landsat 7, 8 and future Landsat 9, a characterisation of temporal dynamics in biophysical parameters becomes possible at decametric resolution. Accurate retrieval and a universal processing workflow for potential Sentinel-2 and Landsat harmonised biophysical products remains a challenge. Previous studies identified MLRAs combined with RTMs in hybrid retrieval workflows as a potential solution. This study investigated the impact of multiple properties of such workflows on the retrieval performance for LAI over a Dutch beech forest site.

Addition of AI noise on the RTM spectral database was found to be most important for prediction performance with a difference of 1.09 m2 m−2 in median RMSE compared to no noise. A level of 5% was optimal in this study. On the other hand, MI noise showed less improvements. Added noise helps the MLRAs to generalise and prevent over-fitting on the pure RTM output. Previous studies did not investigate the effect of different noise definitions and some did not report precisely how noise was defined. With respect to its importance, a clear definition and careful sensitivity analysis should be paramount for future studies.

The choice for the heterogeneous DART RTM in comparison with the turbid medium PROSAIL model resulted in a median RMSE difference of 0.42 m2 m−2. An additional advantage of DART in this study was the lower spread of performance of other inversion workflow features, i.e., DART led to more consistent inversion workflows. Apart from this, care must be taken when using the PROSAIL model with biochemical prior information, as PROSAIL is very sensitive to this. The choice of a specific MLRA was found to be less critical in terms of prediction performance. However, MLRAs varied significantly in run-time, also depending on the implementation and code optimisation. When choosing a particular MLRA, these secondary benefits should be weighted together with the expected accuracy.

**Author Contributions:** Conceptualization, B.B., J.V. (Jochem Verrelst), J.G.P.W.C.; methodology, B.B., J.V. (Jochem Verrelst), J.-P.G.-E., J.G.P.W.C.; formal analysis, B.B., J.-P.G.-E., L.G., J.d.O.; data curation, B.B., L.G.; writing–original draft preparation, B.B., J.G.P.W.C.; writing–review and editing, B.B., J.V. (Jochem Verrelst), J.-P.G.-E., J.G.P.W.C., J.V. (Jan Verbesselt); visualization, B.B.; supervision, J.G.P.W.C., J.V. (Jan Verbesselt), M.H.; project administration, M.H.; funding acquisition, M.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the IDEAS+ and QA4EO contracts funded by ESA-ESRIN. Jochem Verrelst was funded by the European Research Council (ERC) under the ERC-2017-STG SENTIFLEX project (grant agreement 755617).

**Acknowledgments:** Data analysis was supported by the Research and User Support (RUS) Service. The RUS Service is funded by the European Commission, managed by the European Space Agency, and operated by CSSI and its partners. The authors thank the Dutch Forestry Service (Staatsbosbeheer) for granting access to the site, and Devis Tuia and Benjamin Kellenberger for helpful discussions on machine learning. We thank two anonymous reviewers for their valuable comments.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A**

**Figure A1.** Comparison of Fieldspec inverted *Cab* with SPAD meter value taken from the same leaf samples. Solid line, *R*2, *p*-value and formula correspond to the quadratic fit. Dates correspond to the leaf sampling dates.

**Table A1.** List of TLS and satellite sampling dates.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Assessment of UAV-Onboard Multispectral Sensor for Non-Destructive Site-Specific Rapeseed Crop Phenotype Variable at Di**ff**erent Phenological Stages and Resolutions**

#### **Sadeed Hussain 1, Kaixiu Gao 1, Mairaj Din 2, Yongkang Gao 1, Zhihua Shi 1,3 and Shanqin Wang 1,3,\***

Received: 7 January 2020; Accepted: 22 January 2020; Published: 26 January 2020

**Abstract:** Unmanned aerial vehicles (UAVs) equipped with spectral sensors have become useful in the fast and non-destructive assessment of crop growth, endurance and resource dynamics. This study is intended to inspect the capabilities of UAV-onboard multispectral sensors for non-destructive phenotype variables, including leaf area index (LAI), leaf mass per area (LMA) and specific leaf area (SLA) of rapeseed oil at different growth stages. In addition, the raw image data with high ground resolution (20 cm) were resampled to 30, 50 and 100 cm to determine the influence of resolution on the estimation of phenotype variables by using vegetation indices (VIs). Quadratic polynomial regression was applied to the quantitative analysis at different resolutions and growth stages. The coefficient of determination (R2) and root mean square error results indicated the significant accuracy of the LAI estimation, wherein the highest R2 values were attained by RVI = 0.93 and MTVI2 = 0.89 at the elongation stage. The noise equivalent of sensitivity and uncertainty analyses at the different growth stages accounted for the sensitivity of VIs, which revealed the optimal VIs of RVI, MTVI2 and MSAVI in the LAI estimation. LMA and SLA, which showed significant accuracies at (R2 = 0.85, 0.81) and (R2 = 0.85, 0.71), were estimated on the basis of the predicted leaf dry weight and LAI at the elongation and flowering stages, respectively. No significant variations were observed in the measured regression coefficients using different resolution images. Results demonstrated the significant potential of UAV-onboard multispectral sensor and empirical method for the non-destructive retrieval of crop canopy variables.

**Keywords:** unmanned aircraft vehicle; multispectral sensor; vegetation indices; rapeseed crop; site-specific farming

#### **1. Introduction**

Applications of the unmanned aerial vehicle (UAV) platform in precision agriculture (PA) offered a precise and reliable solution for the optimisation of crop monitoring and management [1]. UAV in PA is a new but reliable remote sensing tool that can acquire high spectral, spatial and temporal resolution data and has the advantages of low cost, flexible platform and bird's eye view for rapid data collection [2,3]. UAV in PA can customise the imaging sensor to meet the spectral and spatial requirements and to achieve fast utilisation [4]. Compared with other remote sensing platforms, UAVs can fly at low altitudes and can capture high-resolution imagery, which offers very detailed spectral and spatial

descriptions of the field [1,4]. For highly dynamic vegetation monitoring, the UAV-onboard sensor allows the precise determination of plant location during the growing season to act with corrective measures [2,4,5]. The timely monitoring of crop dynamics is critical in addressing various forms of environmental issues, agricultural practices, degasses and pest control that can eventually lead to enhanced production [6–8]. UAV, sensors and other geospatial techniques (e.g., GIS, remote sensing and GPS) help in determining field variations and optimum fertilisation, disease diagnosis and pest control for sustainable production; this type of determination process is called PA [2,9]. The application of remote sensing in agriculture has become an important research direction because it provides valuable information in terms of agronomic parameters, which facilitate sustainable crop monitoring. Remote sensing data also produce repeated and useful non-destructive crop biophysical attributes for PA [7,10]. Gathering these information required for the enhancement of agricultural production in PA is the key issue that requires the use of a proper design and decision-making system [2,11].

The data collected through the aerial and satellite approaches of remote sensing are useful; however, these platforms are limited by temporal and spatial resolutions, cloud covers and operational costs [12,13]. During the last decade, UAV has been typically used in crop monitoring with high-spatial resolution imagery, low operational cost and near-real-time data acquisition based on powerful sensor-bearing ideal platforms for mapping and monitoring [2,14,15]. In the practical applications of remote sensing techniques in PA, various agricultural equipment units have been developed but require major improvements to meet the variation requirements of seasonal patterns in numerous indices that are used in agronomic evaluations [7].

Leaf area index (LAI), leaf mass per area (LMA) and specific leaf area (SLA) are leaf functional traits that provide information about vegetation canopies in the functional diversity assessment and quantification of physiological processes to prescribe optimum management strategies [16]. These traits correspond to the analytical variables of the plant physiological process, such as growth rates, photosynthetic capacity and plant life strategies [17]. LAI refers to the one-sided green leaf area per unit of ground surface area and is a dimensionless parameter that characterises photosynthesis and respiration and defines the functional link to the canopy spectral reflectance [8,10,18]. LMA is calculated as the ratio of leaf dry mass to leaf area, which is a key biophysical variable involved in plant light capture and carbon gain [19]. Numerous studies identified its significance to the photosynthesis–nitrogen relationship and other essential characteristics, such as the leaf mass-based nitrogen content, which can be inferred from the remote sensing data as the leaf dry mass composed of several organic elements that absorb radiation at a fixed wavelength [19,20]. SLA refers to the leaf area per unit of dry leaf mass or the leaf mass per unit area, which are usually expressed in m2/kg. SLA is an important trait associated to the plant growth rate, which provides information about photosynthetic capacity and leaf nitrogen variations [17]. This trait is indicative of the physiological processes, such as growth rate, light capture and survival [21].

Crop growth and yield models require the estimates of crop biophysical parameters and the association of infrared reflectance to the biophysical parameter of crop canopies that offer the tool for linking multispectral remote sensing to crop growth behaviour [8]. For instance, the soil and water assessment tool model uses the maximum LAI as the key parameter that influence crop growth and watershed flow routing; this parameter can also be changed according to the agronomists' feedback [22,23]. Agricultural studies need to compute and monitor the biochemical and biophysical properties of plants for crop growth, chlorophyll content, nitrogen content and LAI that provide clues to crop health and productivity. Compared with direct field surveys, the remote sensing techniques for the acquisition of these parameters are performed in real time, are non-destructive and provide the spatial details for the measurement and monitoring of these parameters. The empirical relationship between the biophysical parameters and spectral vegetation indices (VIs) has been developed and used to address the relationship between crop traits and canopy reflectance. The empirical relationship of spectral data is effective in predicting the crop biophysical traits because of its simplicity and ability to capture various canopy variations, and this approach is widely used because of its ease in computation [24].

VIs have been developed as a mathematical combination of various spectral bands of the electromagnetic spectrum and are related to various canopy parameters, which can enhance the vegetation signals by reducing the atmospheric and soil effects. The spectral VIs measure vegetation activity and demonstrate the spatial variations of the different seasons through space [7,25]. Substantial developments have been made in the interpretation and analysis of the spectral VIs, which have been applied from the field to the global level [25]. However, these indices are species-specific and are therefore not robust in different species with various leaf and canopy architectures [2,7,25].

Several studies regarding the application of UAV in rice crop yield and biomass estimation have been published [13]. The UAV application in sunflower crop has shown a significant correlation between normalized difference vegetation index (NDVI) and grain yield biomass and applied nitrogen content [2], whereas the UAV platform ('VIPtero') for onboard multispectral camera in vineyard management demonstrated a clear crop heterogeneity, which was consistent with the ground observations [9]. The multispectral sensor mounted on a UAV over potato crops demonstrated a significant correlation between NDVI and green normalized vegetation index (GNDVI) with LAI, plant cover and chlorophyll content [12].

Although the application of UAV has progressed considerably, it still has many limitations, such as payload, cost and operation consistency in agriculture. Thus, a substantial amount of work is still required. UAVs capture numerous high-resolution images that require lengthy processing times, high storage devices and high labour intensity even for professionals because of its low-altitude image acquisition capability [26,27]. Fine-resolution VIs demonstrated a significant influence on the accuracy of prediction models [28]. Although studies have determined that high-resolution images might not be the optimal choice for environmental variable and concluded that fine-resolution VIs obtain a similar prediction model accuracy [2,29,30]. The measuring scale is dependent on the association between biophysical traits and VIs from images with different resolutions, which can help select the appropriate resolution for site-specific management.

Remotely sensed VIs can offer significant information derivatives to agronomic parameters at the field scale. Therefore, the sensor that can monitor crops during the growth season at high spectral, spatial and temporal resolutions will provide useful information in the efficient and sustainable crop management with site-specific basis. Amongst the phenotypical parameters derived using remotely sensed data, LMA and SLA have received minimal attention. However, the retrieval of the mentioned parameters using the radiative transfer model (RTM) inversion could be ill-advised and computationally demanding because it requires a number of leaf and canopy variables. In this study, we optimised and validated numerous VIs for LAI and leaf dry weight (DW) prediction and investigated the potential of empirically estimated LAI and DW for the estimation of the two leaf functional traits, namely, LMA and SLA. Therefore, our study mainly aims to (1) evaluate the potential of a multispectral sensor onboard a UAV in relation to crop biophysical parameters, such as LAI and DW, during the growing season over rapeseed crops at different phenological stages, (2) assess the spectrally predicted leaf DW and LAI for the calculation of LMA and SLA at different phenological stages and (3) assess the effects of resolution on the LAI prediction over the phenological stages to provide simple, rapid and useful information for PA.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The field experiment was conducted at the Zishi experimental station of Huazhong Agricultural University in Jingzhou (30◦11 28.62" N, 112◦23 31.61" E) located in the south-central part of Hubei Province, China from 18 December 2017 to 10 March 2018. This area has a subtropical monsoon climate with an average rainfall of 1100–1300 mm, average annual temperature of 15.9 ◦C to 16.6 ◦C, hydrothermal synchronisation and consistent good climatic condition suitable for agricultural activities. Winter rapeseed was used as the experimental material in this study. To imitate the high variability of growing conditions and investigate the spectral VIs for the estimation of crop phenotype parameters,

the data were collected from 16 different fields and 170 point locations. The experimental fields and sampling points are presented in Figure 1.

**Figure 1.** Zishi experimental station and field view.

#### *2.2. Biophysical Parameter Measurements*

The non-destructive LAI was measured using a plant canopy analyser (SunScan, Probe type SS1, Delta-T Devices, England) over all point locations (Figure 1) after 55 days of plantation on 18 December 2017, after 79 day of plantation on 11 January 2018, 106 days after plantation on 7 February 2018 and after 138 day of plantation on 10 March 2018 according to the crop growth stages from seedling to maturity during the growing season. Measurements were obtained on clear cloudless days. At each phenological stage, 30 subsampling points were selected randomly, and four plants were sampled destructively with their roots. The plant density was also recorded in each point, location, latitude and longitude. The samples were stored in bags and moved to the laboratory. The leaves, stems and roots of the samples were separated to obtain their fresh weights and then placed in an oven at 70 ◦C for 48 hours to obtain the dry mass. The fresh and dry weights of each leaf sample were determined using a high-precision digital scale. LMA and SLA were calculated as follows:

$$LMA = Dw/LA \tag{1}$$

where Dw and LA are the leaf DW and leaf area, respectively.

$$SLA = LA / Dw\tag{2}$$

where LA and Dw are the leaf area and the corresponding Dw, respectively [17,31].

#### *2.3. UAV Image Acquisition*

Multispectral images were also acquired on the days of LAI measurement using the UAV platform. Prior to the flight campaign, the camera was mounted on the UAV under suitable working condition. A lightweight RedEdge MicaSense multispectral camera (MicaSense RedEdge®) onboard a fixed-wing UAV (T-EZ; Golden Wing UAS Co., Ltd.; Chengdu, China) was used to capture the field images shown in Figure S1. A stable platform was used to adjust the camera pointing towards the nadir. The spectral bands covered the wavelength intervals 450 nm centre, 20 nm bandwidth (blue), 560 nm centre, 20 nm bandwidth (green), 668 nm centre, 10 nm bandwidth (red), 717 nm centre, 10 nm bandwidth (RedEdge), 840 nm centre and 40 nm bandwidth (near-infrared). After take-off, the UAV was programmed to follow the predefined route and complete the flight campaign. The flight altitude was 300 m and resulted in a ground resolution of 20 cm.

#### *2.4. Image Processing*

The raw multispectral images acquired by the RedEdge camera were processed in terms of image mosaicking, vignetting correction and raw digital number (DNs) for the reflectance conversion using the Pix4Dmapper (Pix4D SA, Switzerland) software. The raw digital numbers of the acquired images were firstly converted into radiance and then into surface reflectance based on the linear regression models by using the calibration target data, wherein the surface reflectance is a linear function of the digital numbers [32]. In each stage, 1024 images were captured and mosaicked from the study area, and the georeferenced images of the study area were acquired at the end of the fourth crop season. Moreover, the raw images with the ground resolution of 20 cm were resampled to 30, 50 and 100 cm by using the pixel aggregate interpolation method to assess the influence of image resolution on VIs. Lastly, the VIs were calculated.

#### *2.5. Multispectral VIs*

The canopy spectral data were used to develop the VIs, which are sensitive to canopy structure, pigments and chlorophyll, for the LAI, LMA and SLA estimation in rapeseed crops. Nine VIs were calculated, several of which have been proposed as surrogates for LAI estimation [10,25,32]. The calculated VIs, their formula and description are listed in Table 1. For all VIs, the images with the same resolution were stacked, and the mean VI values centred on the 1 m2 location were extracted for 160 point locations in ENVI 5.1 (Exelis Visual Information Solution, Inc.; Boulder, CO, USA).

**Table 1.** Description and formulas of the investigated VIs.

Note: RVI, ratio vegetation index; NDVI, normalized difference vegetation index; GNDVI, green normalized vegetation index; BNDVI, blue normalized difference vegetation index; SAVI, soil-adjusted vegetation index; OSAVI, optimised soil-adjusted vegetation index; MSAVI, modified soil-adjusted vegetation index; MSAVI2, modified soil-adjusted vegetation index 2; MTVI2, modified triangular vegetation index 2.

#### *2.6. Statistical Analysis*

Regression models were used for the quantitative analysis between spectral VIs as independent and LAI and leaf DW as dependent variables at the different phenological stages. A total of 140 and 30 sample points for LAI and LD were used for model calibration, respectively, whereas 30 sample points each were used for validation. To assess the model performance, the coefficient of determination (R2), root mean square error (RMSE) and relative RMSE (RRMSE) were used and calculated using Equations. (3), (4) and (5), respectively [18,42,43]. The large R2 value and low RMSE and RRMSE values indicated the high precision and accuracy of the model as follows:

$$R^2 = 1 - \frac{\sum\_{i=1}^{n} \left( y\_i - \hat{y}\_i \right)^2}{\sum\_{i=1}^{n} \left( y\_i - \overline{y} \right)^2} \tag{3}$$

$$RMSE = \sqrt{\frac{\sum\_{i=1}^{n} \left(y\_i - \hat{y}\_i\right)^2}{n}} \tag{4}$$

$$\text{RRMSE} = \frac{\text{RMSE}}{\overline{y}} \times 100\tag{5}$$

where R2 and RMSE measure the relationship of the fitting function. However, for the nonlinear functions amongst LAI and VIs, the R2 values may be misleading because sensitivity changes continuously [10]. Thus, a different accuracy matric noise equivalent (NE) was needed to evaluate and determine the LAI estimation accuracy. NE also accounts for the slope and scattering of data points in the best-fit function as follows [10,25,44]:

$$\text{ME}\DeltaLAI = \frac{RMSE \{VI \, VS.LAI\}}{d(VI)/d(LAI)}\tag{6}$$

where RMSE is the root mean square error of the best-fit function between the spectral Vis and LAI, and d(VI)/d(LAI) is the first-order derivative of the observed relationship.

#### **3. Results**

#### *3.1. Association of VIs to the Phenolgical Stages*

A wide range of variations in LAI, DW, LMA and SLA due to different canopy architectures as well as variations due to leaf sizes and shapes at different phenological stages was recorded. Table 2 presents the descriptive statistics of canopy LAI, DW, LMA and SLA of rapeseed oil at the phenological stages in the Zishi experimental station.

All the VIs investigated at the different growth stages revealed close associations with the crop reflectance characteristics. The spectral VIs demonstrated significant variations in all the phenological stages from seedling to maturity, and the results are compiled in Table 3. The VIs attained its peak values at the elongation stage and started to decline afterwards. All the VIs demonstrated a significant relation with LAI at the seedling stage. Furthermore, the MSAVI2, MTVI2 and SAVI showed the most sensitive VIs with an R2 value of 0.87, followed by MSAVI (R2 = 0.86), whilst NDVI and OSAVI achieved R2 = 0.85, GNDVI R2 = 0.82 and BNDVI obtained an R2 value of 0.81 at the seedling stage. RVI, NDVI, MSAVI and MTVI2 showed the most influenced VIs to LAI from the seedling to maturity stages and attained the maximum R2 (0.93–0.89, 0.88–0.86, 0.87–0.86 and 0.89–0.86) from the elongation to flowering stages. The minimum R2 values at the maturity stage were expected due to crop senescence.

**Table 3.** Relationship of spectral VIs at the different phenological stages.

#### *3.2. Evaluation of Spectral VIs for Estimation of Rapeseed LAI*

The optimal spectral VIs for the estimation of the LAI best-fit model relationship for VIs and LAI at all the phenological stages are shown in Table 4. The spectral pattern variations and close association with crop reflectance revealed the potential capability of the explored VIs in LAI estimation. All the VIs demonstrated a high significant relationship to the temporal distribution of LAI at the different growth stages, except at the maturity stage, which revealed the lowest coefficient of determination (R2 < 0.50) and were omitted from the calculations. Amongst the investigated phenological stages, all the VIs attained their highest coefficient of determination at the elongation stage, whereas the highest R2 value was achieved by RVI (R2 = 0.93), followed by MTVI2 (R2 = 0.89) at the different phenological stages, as provided in Table 3. At the next flowering phenological stage, R2 started to decrease until maturity, and the best-fit model of LAI with spectral VIs provided RVI (R2 = 0.89), followed by NDVI (R2 = 0.87) and MTVI2, MSAVI (R2 = 0.86) and OSAVI (R2 = 0.84), whereas SAVI attained R2 = 0.83 and GNDVI and MSAVI2 attained R2 = 0.82. The least R2 value was attained by BNDVI (R2 = 0.79). MSAVI and MTVI2 have the modifying factor to deal with the saturation problem when the LAI surpasses the saturation level (LAI > 3), whilst the sensitivity of NDVI suffered and levelled out.

**Table 4.** Quadratic polynomial regression and summary statistics of relationship between LAI and studied VIs based on calibration data set at different phenological stages.

#### *3.3. Sensitivity Analysis*

The sensitivity analysis results demonstrated that BNDVI and GNDVI presented the maximum insensitivity, whereas OSAVI exhibited moderately higher sensitivity and was unreliable during the entire range of LAI. NDVI and SAVI exhibited higher sensitivity for LAI < 2.5 m2m−2 and a decreasing sensitivity trend was observed for higher LAI > 3.0 m2m−2, whereas MSAVI2 performed better than NDVI and SAVI even at LAI > 3.0 m2m−2 [25,45]. RVI showed higher sensitivity to LAI and is the best index for the detection of numerical variations in LAI because other indices are not robust at higher LAI values. MSAVI and MTVI2 exhibited higher and consistent sensitivity for LAI < 3.0 m2m−2, and a decreasing trend was observed for higher LAI > 3.0 m2m−2 for MSAVI. Therefore, RVI and MTVI2 demonstrated the lowest NE values to ensure the highest sensitivities to LAI and the best multispectral indices for quantitatively detecting variations in LAI (Figure 2).

**Figure 2.** Sensitivity analysis for VIs evaluated for LAI estimation.

#### *3.4. Relationship of VIs and Leaf DW*

To evaluate the optimal VIs in the estimation of leaf DW (gm−2), the best-fit model relationship was used between the leaf DW and VIs presented in Table 1 at the different phenological stages demonstrated in Table 5. The best-fit function demonstrated the nonlinear relationship between VIs and leaf DW with R2 that range from 0.34 to 0.68 and maximum R2 (0.68) attained by using SAVI at the seedling stage. At the early seedling stage, SAVI performed well with the maximum R2 value, which was demonstrated at the early growth stage, wherein the canopy is not fully developed, and the soil contributed to the reflectance. The highest R2 (0.75) value was obtained by RVI and GNDVI, followed by NDVI with R2 = 0.70 at the elongation stage and then began to decline afterwards. At the flowering stage, the maximum R2 (0.72) value was achieved by MSAVI, followed by NDVI, SAVI and MSAVI2 at 0.71. All the VIs showed similar asymptotic patterns with a certain degree of scattering from the nonlinear fits.

**Table 5.** Power function models and summary statistics of the relationship between LD and the explored VIs based on the calibration data set at different phenological stages.

#### *3.5. Estimation LMA and SLA using Spectral VIs*

The performance of all the indices were compared at the different phenological stages to select the optimum index in the prediction of LAI and DW. The best-fit regression model illustrated the good relationship in the prediction of LAI and DW (Tables 4 and 5), respectively, at nearly all growth stages, except the maturity stage, which was excluded from the data. The best-fit model equations of the optimal VIs were applied for the prediction of LAI and DW at each phenological stage and used in Equations (1) and (2) in the estimation of LMA and SLA, respectively. The spatial distributions of LAI, LMA and SLA are shown in Figures S2–S4. The optimal prediction equation was used in the LAI estimation, whereas LMA and SLA were utilised in Equations (1) and (2) for the generation of maps at all the explored phenological stages, respectively. Field measurement with the estimations from composited maps provided the validation results and illustrated the potential use of LAI, LMA and SLA derivation at all the growth stages.

#### *3.6. Validation of LAI, LMA and SLA Estimates*

To validate the calibration model for LAI and estimated LMA, SLA with the predicted LAI and DW, multispectral images and ground measured data were used whilst the determined optimal VIs were considered in the prediction of LAI and DW at each growth stage and used in the estimation of LMA and SLA. The observed and estimated LAI, LMA and SLA were compared in Figures 3–6. The predictions of optimal VIs were selected on the basis of the highest R2 and lowest RMSE, which were plotted against the ground observed LAI, LMA and SLA. The optimal VIs, RVI and MTVI2 suggested consistency between the measured and estimated LAI (R2 > 0.80) of rapeseed oil crops at all the growth stages. By contrast, the model validation demonstrated less adequate results for LMA and SLA at the seedling stage (R2 = 0.62, 0.25) but realised consistent results at the elongation (R2 = 0.85, 0.80) and flowering (R2 = 0.85, 0.71) stages (Figure 6). Plotting the measured–estimated LMA values demonstrated less deviation from the 1:1 line, whereas SLA was strong under and over the prediction that occurs at higher SLA values, thereby suggesting the saturation of the spectral signal [46].

**Figure 3.** Estimated vs. measured LAI at the seedling stage. The green line shows the 1:1 correlation of the estimated and measured variables, whereas the red lines present the linear regression models.

**Figure 4.** Estimated vs. measured LAI at the elongation stage. The green line shows the 1:1 correlation of the estimated and measured variables, whereas the red lines present the linear regression models.

**Figure 5.** Estimated vs. measured LAI at the flowering stage. The green line shows the 1:1 correlation of the estimated and measured variables, whereas the red lines represent the linear regression models.

**Figure 6.** Estimated vs. measured LMA and SLA at the (**a**) seedling, (**b**) elongation and (**c**) flowering stages. The green line shows the 1:1 correlation of the estimated and measured variables, whereas the red lines are the linear regression models.

#### *3.7. Evaluation of Image Resolution E*ff*ect*

Four images with 20 cm pixel resolution were acquired over the rapeseed crops during the growing season and resampled to 30, 50 and 100 cm pixel resolutions. To evaluate the effects of image resolutions on the remote estimation of LAI over the rapeseed crops, the quantitative relationship between VIs and LAI were obtained at different image resolutions (20, 30, 50 and 100 cm). At all the mentioned image resolutions, the estimated LAI using VIs had no significant differences. The precision of the best-fit regression model between the LAI and VIs determined no significant variations and was similar at all the image resolutions. Thus, a significant high correlation between spectral VIs and LAI was acquired at the elongation and flowering stages, which are independent of image resolution. The best-fit model R2 and RMSE are presented in Figure 7.

**Figure 7.** Best-fit model R2 and RMSE that describe the relationship between spectral VIs and LAI with respect to the image resolution at the (**a**) seedling, (**b**) elongation, (**c**) flowering and (**d**) maturity stages.

#### **4. Discussion**

As a remote sensing platform, UAVs are ideal for mapping and monitoring in PA [9,15]. With the recent developments in UAV technologies, numerous studies have been conducted using remote sensors onboard UAVs to evaluate its application for PA and subsequently obtained considerable attention from researchers [2,9,12–15,25,31,47].

The key canopy biophysical and biochemical parameter retrievals can be attained adequately using certain spectral bands, that is, the green spectrum is ideal for leaf chlorophyll estimation [35]. Reflectance in the near-infrared (NIR) is recognised as sensitive to variant biomass and canopy structures, whereas the red and green reflectance values respond to the background variations and senescent vegetation [48].

This study aimed to determine the optimal VIs derived from UAV multispectral images that could be used for the remote and empirical non-destructive estimation of biophysical parameters, such as LAI, LMA and SLA over rapeseed crops at the different growth stages. The empirical method for the estimation of biophysical parameters is proficient and delivers precise estimations [24,32]. Moreover, this type of method has been widely applied for the improved quantitative accuracy of spectral VIs and used for the evaluation of crop attributes, such as LAI, N status and biomass for various crops [10,49]. In this study, the potential capability of multispectral VIs mentioned in Table 1 was evaluated over rapeseed crops at different growth stages. The results revealed that the crop variations from seedling to maturity are similar to the canopy reflectance characteristics. Variations in the crop canopy reflectance attains its peak value at the elongation stage and begins to decline afterwards until the maturity stage; this behaviour is due to crop senescence, as reported previously in [10,50]. Variations in canopy reflectance at the different spectral wave bands are synchronised with LAI because of the modification of leaf chlorophyll contents at different phenological stages [51]. Spectral reflectance from the canopy is affected not only by the biophysical characteristics of foliage but also by the direction of incidence radiation, canopy architecture and soil background [52]. The small variations in spectral reflectance at the early phenological stage in the visible region are likely due to the crop nitrogen content and soil water background. However, variations in NIR reflectance with variant leaf orientations at certain growth stages are due to the overlapping leaves that decrease the active photosynthetic size as the LAI reaches a plateau [53,54].

Table 2 demonstrates the highly significant relationship between LAI with all the aforementioned spectral VIs derived from the multispectral remote sensing data. However, a close relationship was obtained by RVI with LAI because LAI had a value of approximately 3 m2m−2, followed by NDVI and MTVI2 [55]. When LAI reached its saturation point, the reduced variability of red and NIR reflectance make RVI and NDVI insensitive. However, GNDVI apparently did not reach the saturation level even at LAI with moderate to high (4–5) values. The precision of the LAI estimate was better when green and blue bands were used instead of that that with the red band even at LAI values greater than 3 m2m−2 [42]. The NIR band has a considerably strong impact on the relationship between LAI and spectral reflectance and must be considered under varying crop situations [7,10,56]. The nonlinear best-fit function between LAI and VIs demonstrated inconsistent sensitivity, and R2 and RMSE values for estimating the accuracy of LAI can be misleading. Thus, the sensitivity analysis (NE) can be applied as the precision indicator to verify the performance of VIs in LAI estimation [57,58].

LMA and SLA are important in determining plant composition and eco-physiological characterisation [17,59–62]. Compared with the general empirical approach for the estimation of SLA, as reported in [32,48], its indirect estimation using the predicted LAI and DW illustrated highly significant retrieval accuracy. The SLA results showed remarkable variations between the different growth stages, as shown in Figure 6. The estimation results of SLA at the elongation and flowering stages are consistent with those of previously reported studies [63]. Soil background has less affected the correlation results as the data were captured at mature seedling (rosette) stage. The slight change in the performance of multispectral VIs from the UAV images could be awarded to the different viewing geometries and instrument spectral response functions [32].

Remotely sensed LMA (R2 = 0.85, nRMSE= 0.10) and (R2 = 0.85, nRMSE = 0.12) at the elongation and flowering stages was calculated with high estimation accuracy, as previously reported in [47]. The study estimated LMA using the LUT-based PROSAIL inversion method over the rapeseed crop. The LMA estimation results were also consistent with [19] who estimated LMA across a wide range of plant species through continuous wavelet analysis, and the study demonstrated remarkable implications in the prediction of dry matter content at the canopy level.

In terms of image resolution, our study concluded that the application of very fine-resolution remote sensing images does not reflect any significant difference. Thus, the fine-resolution images would not essentially increase the prediction accuracy of the regression model between LAI and spectral VIs. The precision of regression model between LAI and VIs did not change significantly by using images with different resolutions, as shown in Figure 7. The end users of UAV images need to consider the spatial, spectral and temporal resolutions and processing time for site-specific crop monitoring and management, as previously reported in [2,29,64–66].

#### **5. Conclusions**

In summary, herein we reported the application of multispectral sensors onboard UAV in acquiring images of rapeseed crops and retrieving their biophysical characteristics, such as LAI, LMA and SLA, throughout the phenological stages. The results revealed that the strongest relationship of spectral VIs was obtained in the elongation stage. The sensitivity analysis between LAI and VIs revealed that RVI, MSAVI and MTVI2 were the optimal VIs for LAI estimation. LMA and SLA results demonstrated the significant estimation accuracy by using the predicted leaf DW and LAI at the elongation and flowering stages, respectively. In terms of image resolutions, robust results can still be obtained when the maximum 100 cm resolution imagery is used for oil rapeseed crop characterisation. Therefore, high-altitude images will be preferred in obtaining a decreased number of images, which will significantly influence image acquisition and processing time.

Compared with earlier remote sensing platforms, the sensor onboard UAVs had the basic advantage of reaching the targeted site and acquiring very comprehensive information throughout the growing season. Unlike satellite- and aircraft-based remote sensing platforms, the UAV remote sensing platform demonstrated better operational advantages, such as the provision of high-spatial resolution imagery with low operational cost and near-real-time data acquisition. Although satellite data were limited because of temporal resolution, cloud cover, availability at the ideal time and the operational cost for high-resolution imagery, this technology has become an important tool for large area mapping and monitoring because of its distinctive capabilities in satellite remote sensing imagery that offers extensive information in a synoptic and frequent manner.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2072-4292/12/3/397/s1, Figure S1: The UAV used for image acquisition, Figure S2: Map of LAI, LMA and SLA at seedling stage, Figure S3: Map of LAI, LMA and SLA at elongation stage, Figure S4: Map of LAI, LMA and SLA at flowering stage.

**Author Contributions:** Conceptualization, S.S. and S.W.; Data curation, S.S.; Formal analysis, S.S., K.G., M.D. and Y.G.; Funding acquisition, S.W.; Investigation, K.G. and S.W.; Methodology, S.S.; Resources, Z.S. and S.W.; Software, S.H.; Supervision, Z.S. and S.W.; Validation, S.H.; Writing—original draft, S.H.; Writing—review and editing, Z.S. and S.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Key Research and Development Program of China (Grant No. 2018YFD0200900).

**Acknowledgments:** We would like to thanks and acknowledge management of Zishi experiment station, Jingzhou, China and special thanks to Muhammad Yaseen, Muhammad Shakir, Muhammad Asim, Nauman and Sajjad for their help at various stages of experiments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Turner, D.; Lucieer, A.; Watson, C. An Automated Technique for Generating Georectified Mosaics from Ultra-High Resolution Unmanned Aerial Vehicle (UAV) Imagery, Based on Structure from Motion (SfM) Point Clouds. *Remote Sens.* **2012**, *4*, 1392–1410. [CrossRef]

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Retrieval of the Fraction of Radiation Absorbed by Photosynthetic Components (***FAPARgreen***) for Forest Using a Triple-Source Leaf-Wood-Soil Layer Approach**

**Siyuan Chen 1,2, Liangyun Liu 1,\*, Xiao Zhang 1,3, Xinjie Liu 1, Xidong Chen 1,3, Xiaojin Qian 1,3, Yue Xu 4 and Donghui Xie 4**

Received: 12 September 2019; Accepted: 18 October 2019; Published: 23 October 2019

**Abstract:** The fraction of absorbed photosynthetically active radiation (FAPAR) is generally divided into the fraction of radiation absorbed by the photosynthetic components (*FAPARgreen*) and the fraction of radiation absorbed by the non-photosynthetic components (*FAPARwoody*) of the vegetation. However, most global FAPAR datasets do not take account of the woody components when considering the canopy radiation transfer. The objective of this study was to develop a generic algorithm for partitioning *FAPARcanopy* into *FAPARgreen* and *FAPARwoody* based on a triple-source leaf-wood-soil layer (TriLay) approach. The LargE-Scale remote sensing data and image simulation framework (LESS) model was used to validate the TriLay approach. The results showed that the TriLay *FAPARgreen* had higher retrieval accuracy, as well as a significantly lower bias (R2 = 0.937, Root Mean Square Error (RMSE) = 0.064, and bias = −6.02% for black-sky conditions; R2 = 0.997, RMSE = 0.025 and bias = −4.04% for white-sky conditions) compared to the traditional linear method (R2 = 0.979, RMSE = 0.114, and bias = −18.04% for black-sky conditions; R2 = 0.996, RMSE = 0.106 and bias = −16.93% for white-sky conditions). For FAPAR that did not take account of woody components (*FAPARnoWAI*), the corresponding results were R2 = 0.920, RMSE = 0.071, and bias = −7.14% for black-sky conditions, and R2 = 0.999, RMSE = 0.043, and bias = −6.41% for white-sky conditions. Finally, the dynamic *FAPARgreen*, *FAPARwoody*, *FAPARcanopy* and *FAPARnoWAI* products for a North America region were generated at a resolution of 500 m for every eight days in 2017. A comparison of the results for *FAPARgreen* against those for *FAPARnoWAI* and *FAPARcanopy* showed that the discrepancy between *FAPARgreen* and other FAPAR products for forest vegetation types could not be ignored. For deciduous needleleaf forest, in particular, the black-sky *FAPARgreen* was found to contribute only about 23.86% and 35.75% of *FAPARcanopy* at the beginning and end of the year (from January to March and October to December, JFM and OND), and 75.02% at the peak growth stage (from July to September, JAS); the black-sky *FAPARnoWAI* was found to be overestimated by 38.30% and 28.46% during the early (JFM) and late (OND) part of the year, respectively. Therefore, the TriLay approach performed well in separating *FAPARgreen* from *FAPARcanopy*, which is of great importance for a better understanding of the energy exchange within the canopy.

**Keywords:** the fraction of radiation absorbed by photosynthetic components (*FAPARgreen*); triple-source; leaf area index (LAI); woody area index (WAI); clumping index (CI); Moderate Resolution Imaging Spectroradiometer (MODIS); soil albedo

#### **1. Introduction**

The fraction of absorbed photosynthetically active radiation (FAPAR) is a significant biochemical and physiological variable used in tracing the exchanges of energy, mass, and momentum, and is also widely used in many climate, ecological, biogeochemical, agricultural, and hydrology models [1,2]. FAPAR is, therefore, an important input parameter and widely used in satellite-based Production Efficiency Models (PEMs) [3–6] to estimate gross primary productivity (GPP) or net primary production (NPP).

In general, the FAPAR inversion algorithms could be divided into two types: empirical statistical models based on vegetation indexes and physical methods based on the canopy radiation transfer model. Although the empirical statistical model based on vegetation indexes is relatively simple, involves only a few parameters and has high computational efficiency, it is subject to many uncertainties due to factors such as the atmospheric environment, vegetation type, and quality of remote sensing data. The physically based methods could be further divided into two categories. The first type is the direct inversion method, which uses the canopy radiation transfer model to link FAPAR with the canopy spectra [7–9]. For example, the Moderate Resolution Imaging Spectroradiometer (MODIS) algorithm uses the three-dimensional radiation transmission model to invert FAPAR from the bi-directional reflectance [10–12]. The Joint Research Centre (JRC) FAPAR algorithm is also based on a physical model that uses a continuous vegetation canopy model [13] to link land surface reflectance with FAPAR. However, these methods are mostly based on the radiative transfer model; thus, the inversion process is complicated for retrieval of FAPAR. The main problem with such methods is that it is difficult to overcome the uncertainty caused by model coupling and spatial heterogeneity. The second type of physically based method is the forward modeling method [14–19]. Most models of this type are based on the gap fraction model, which determines FAPAR according to canopy structure parameters such as LAI and the clumping index. The disadvantage of this approach is that it relies too much on the accuracy of the canopy structure parameters. Furthermore, it is difficult to accurately determine the soil albedo and extinction coefficient, which are also important parameters needed to determine the contribution of multiple scattering between the soil and canopy to FAPAR [17].

Recently, several global FAPAR products have become available, including theModerate Resolution Imaging Spectroradiometer (MODIS) [20,21], Energy Balance Residual (EBR) [15], Multi-angle Imaging SpectroRadiometer (MISR) [11], CYCLOPES [22], GLOBCARBON [23], Global Land Surface Satellite (GLASS) [14], the Medium Resolution Imaging Spectrometer (MERIS) [24], Joint Research Center Two Stream Inversion Package (JRC-TIP) [25], and European Space Agency (ESA) products [26]. These global FAPAR products have been widely validated, with reported errors varying from 0.08 to 0.23 [14,27–32]. However, most of the global FAPAR products do not consider the effect of non-photosynthetic components in the radiative transfer process, which introduces errors, especially for forest types. Moreover, many researchers have used the fraction of radiation absorbed by photosynthetic components (*FAPARgreen*) instead of the fraction of radiation absorbed by the canopy (*FAPARcanopy*) to monitor and estimate the light use efficiency (LUE), radiation use efficiency (RUE), and productivity at different temporal scales [33–36].

The forest vegetation ecosystem plays an important role in the global ecosystem. However, quantifying the temporal variation in *FAPARgreen* for a forest ecosystem represents an important challenge for remote sensing and ecology researchers as it is extremely difficult to measure *FAPARgreen* at large scales over plant growing seasons directly. Also, previous studies have shown that the contribution of woody components is relatively large: for instance, Asner et al. [37] found that stems increased *FAPARcanopy* by 10–40%. Therefore, the partitioning of absorbed radiation into photosynthetic and non-photosynthetic parts is very important for better modeling of vegetation photosynthesis and energy exchange within the canopy.

Already, some studies have looked at the estimation of *FAPARgreen* from remote sensing data. Hall et al. [38] estimated *FAPARgreen* using a simple linear relationship between *FAPARcanopy* and *LAIgreen*/*LAItotal* (*LAItotal* denotes the total leaf area index including green and senescent leaves, while *LAIgreen* represents the green leaf area index). However, this simple partitioning is problematic because

the green and woody components within the canopy do not constitute a simple linear mix in terms of radiation transfer. Zhang et al. [39] first retrieved the biophysical and biochemical variables using the modified PROSPECT model coupled with the SAIL-2 model (hereafter called PROSAIL-2 model), and then calculated *FAPARgreen* and *FAPARcanopy* using the forward simulation approach. However, *FAPARgreen* retrieval using the PROSAIL-2 model is relatively complex and needs several physiological and biochemical parameters as model inputs. Gitelson et al. [40] also separated *FAPARcanopy* into photosynthetically active green components (*FAPARgreen*) and non-photosynthetic active components using the ratio *LAIgreen*/*LAItotal* for maize and soybeans. The relationship between vegetation indices and *FAPARgreen* was also used to retrieve *FAPARgreen* [41,42]. Nevertheless, to date, the current *FAPARgreen* products do not take into account the effect of non-photosynthetic components on canopy radiative transfer.

In this study, we aim to develop an operational algorithm for partitioning *FAPARcanopy* into *FAPARgreen* and *FAPARwoody* for forest types. A simple triple-source leaf–wood–soil layer model (TriLay) that describes the radiation transfer within the canopy-soil system is presented. *FAPARcanopy* is first separated into the fraction of PAR absorbed by the canopy for downwelling radiation (*FAPARcanopy*↓) and the fraction of PAR absorbed by the canopy for the upwelling radiation reflected by the soil background (*FAPARcanopy*↑). Then, *FAPARcanopy*↓ and *FAPARcanopy*↑ are further split into the fraction of radiation absorbed by photosynthetic components (*FAPARgreen*) and that absorbed by non-photosynthetic components (*FAPARwoody*) using the TriLay model. Finally, the *FAPARgreen*, *FAPARwoody*, and *FAPARcanopy* products are generated using the MODIS albedo (MCD43A3), LAI (MCD15A2H), land cover (MCD12Q1), clumping index (CI), and soil albedo products based on the TriLay approach, and the discrepancies between different FAPAR products are used to investigate the contributions of woody components to the canopy-absorbed radiation. The partitioning of absorbed radiation into green and woody parts using the TriLay model is done not just to provide *FAPARgreen* and *FAPARwoody*—which is of great importance for better understanding the energy exchange within the canopy. The consideration of woody components should also improve the accuracy of *FAPARgreen* estimates, which is important for better modeling of vegetation photosynthesis.

#### **2. Materials and Methods**

#### *2.1. Satellite Datasets*

In order to produce *FAPARgreen* products, several satellite datasets were utilized in this study, including MODIS LAI products, land cover products, CI products, and soil albedo products from 2017. Data simulated by the LESS model [43] was used to validate the retrieved *FAPARgreen* products. A mid-latitude region (tile h10v05, covering 30.00◦ N–40.00◦ N and 80.00◦ W–104.43◦ W) was selected to investigate the discrepancy between *FAPARgreen* and other FAPAR products because there were abundant forest vegetation types within this MODIS tile.

#### 2.1.1. MODIS LAI/FAPAR Products (MCD15A2H)

MCD15A2H V006 is a MODIS eight-day composite LAI/FAPAR product that includes FAPAR, LAI, and quality control (QC) data with a resolution of 500 m [44]. The main retrieval algorithm for LAI and FAPAR contains a Look-up-Table (LUT) based on a 3D radiation transfer model [21] that uses the atmospherically corrected Red and near-infrared (NIR) Bidirectional Reflectance Function (BRF) [45]. A back-up algorithm based on the empirical relationships between the Normalized Difference Vegetation Index (NDVI) and LAI and FAPAR at the canopy scale is used at the same time. Also, for the biome types and typical conditions that are considered, observed, and modeled spectral BRFs and soil patterns are compared for each pixel. The LAI and FAPAR values that lie within a fixed level of uncertainty are then taken to be acceptable. Finally, averaged values of LAI and FAPAR are used as the eventually retrieved values [20].

#### 2.1.2. MODIS Land Cover Product (MCD12Q1)

The MODIS Land Cover Type Product (MCD12Q1) provides land cover maps with a temporal resolution of one year and a spatial resolution of 500m at a global scale from 2001 until the present; it includes several classification schemes (International Geosphere-Biosphere Programme (IGBP), University of Maryland (UMD), LAI, BIOME-Biogeochemical Cycles (DBC), Plant Functional Types (PFT), FAO-Land Cover Classification System land cover (LCCS1), etc.). The main algorithm used in MCD12Q1 is a supervised classification method (decision tree) combined with a boosting technique [20] based on MODIS reflectance data [46,47]. The International Geosphere-Biosphere Program (IGBP) classification scheme, which defines 17 land cover types, was used in this study.

#### 2.1.3. Global Clumping Index (CI) Product

The clumping index (CI) signifies the characteristics of groups of foliage in the canopy, and thus, it is a significant parameter describing the structure of the vegetation canopy [48]. According to Jiao et al. [49], a method based on the MODIS Bidirectional Reflectance Distribution Function (BRDF) and a linear relationship between the CI and the normalized difference between the angular indexes of the hotspot and dark spot (NDHD) could be utilized to generate CIs within a valid range (0.33 to 1.00); a back-up algorithm is also used to substitute values for the invalid CIs [49]. A global CI dataset with an eight-day temporal resolution and 500 m spatial resolution covering the period from 2002 to the present has been produced by Jiao et al. [49].

#### 2.1.4. Global Soil Albedo Product

A non-linear spectral mixture model (NSM) model proposed by Liu & Zhang (2018) [15,50] was used to retrieve the global visible (VIS) soil albedo. The main idea in the NSM model is the dual-source vegetation–soil layer approach. In this approach, it is assumed that the leaves are located in the upper canopy, while the soil components are found in the lower canopy. Based on this assumption, the canopy albedo can be approximated as a non-linear mixture of the "pure" vegetation and soil parts [15]. For pixels with abnormal values (smaller than 0.02 or greater than 0.3), prior values acquired by a global database of land surface parameters at 1 km resolution (ECOCLIMAP) [51,52] and the yearly composite value were used instead. Finally, the estimated global VIS soil albedo based on the NSM model is obtained, having a spatial resolution of 500 m and a temporal resolution of eight days.

#### *2.2. Data Simulated by the LESS Model*

To quantitatively evaluate the performance of the proposed TriLay approach for estimating *FAPARgreen*, the LargE-Scale remote sensing data and image simulation (LESS) framework model [43] was employed to generate a simulated dataset that covered most of the conditions found in forests. LESS is a ray-tracing based 3D radiative transfer model which can simulate remote sensing data and images over large-scale and realistic 3D scenes (http://lessrt.org/). LESS employs a weighted forward photon tracing (FPT) method to simulate multispectral bidirectional reflectance factor (BRF) or flux-related data (e.g., downwelling radiation) and a backward path tracing (BPT) method to generate sensor images (e.g., fisheye images) or large-scale (e.g., 1 km2) spectral images. The accuracy of LESS is evaluated with other models as well as field measurements in terms of directional BRFs and pixel-wise simulated image comparisons, which shows very good agreement [43,53].

The modeling area is located in the Genhe Forestry Reserve (Genhe) (120◦12 to 122◦55 E, 50◦20 to 52◦30 N), Greater Khingan of Inner Mongolia, Northeastern China. It has a hilly terrain with 75% forest cover, which is mainly composed of Dahurian Larch (Larix gmelinii) and White Birch (Betula platyphylla Suk.). A pure plot of Dahurian Larch (L9) is established and is selected for modeling, with its location in Figure 1. The position, crown width, breast diameter, tree height, and transmittance for trees were entered into the LESS model. A total of eight scenes were constructed, with each individual tree in a scene consisting of branches and leaves. The branches and size of the woody area were kept the same in all the scenes, but the leaf area was changed to produce scenes with different values of the LAI. For each scene, different FAPAR values (including *FAPARcanopy*, *FAPARgreen* and *FAPARwoody*) were calculated under black-sky (the ratio of diffuse light is zero, and nine different solar zenith angles of 0◦–80◦ at 10◦ intervals were set) and white-sky conditions (the ratio of diffuse light is 1). The main input parameters used in the LargE-Scale remote sensing data and image simulation (LESS) model are listed as Table 1. Finally, a total of 72 black-sky simulations and eight white-sky simulations were achieved for the different conditions giving.

**Figure 1.** Locations of the selected Larch plots (L9) in Genhe Forestry Reserve.

**Table 1.** The main input parameters used in the LargE-Scale remote sensing data and image simulation (LESS) model simulations.

#### *2.3. Algorithms for Estimating Global FAPARgreen and FAPARwoody Datasets*

To split the fraction of radiation absorbed by photosynthetic components (*FAPARgreen*) from *FAPARcanopy*, a novel triple-source leaf–wood–soil layer model was proposed for generating global *FAPARgreen* products for forest vegetation types. Figure 2 is a flowchart of the process used to retrieve global *FAPARgreen* products.

First, *FAPARcanopy* is split into two parts: the fraction of PAR absorbed by the canopy for the downwelling radiation (*FAPARcanopy*↓) and that absorbed by the canopy for the upwelling radiation reflected by the soil background (*FAPARcanopy*↑). Then, *FAPARcanopy*↓ and *FAPARcanopy*↑ are further split into the fraction of radiation absorbed by photosynthetic components (*FAPARgreen*) and that absorbed by non-photosynthetic components (e.g., branches and stems, hereafter called *FAPARwoody*). Finally, *FAPARgreen* and *FAPARwoody* can be calculated separately using the TriLay approach.

**Figure 2.** Flowchart illustrating the Triple-source leaf–wood–soil layer (TriLay) method for estimating the fraction of radiation absorbed by photosynthetic components (*FAPARgreen*) and the fraction of radiation absorbed by woody components (*FAPARwoody*).

2.3.1. The Triple-Source Leaf–Wood–Soil Layer Model

A triple-source leaf–wood–soil layer model (TriLay) was developed to model the radiation transfer within the vegetation–soil system—this model is illustrated as Figure 3.

**Figure 3.** Illustration of the triple-source leaf–wood–soil layer model.

In this study, we used the layer approach to illustrate the distribution of leaves and soil in the whole canopy, which is consistent with the approach used in our previous study [15]. The layer approach assumes that the canopy consists only of green components and woody components and that all leaves are found above the canopy. The NSM model [15] was then used to simulate the canopy albedo.

The fraction of PAR absorbed by the canopy, *FAPARcanopy*, can be separated into the fraction of PAR absorbed by the canopy for the downwelling radiation (*FAPARcanopy*↓) and that for the upwelling radiation reflected from the soil background (*FAPARcanopy*↑):

$$\text{FAPAR}\_{\text{canopy}} = \text{FAPAR}\_{\text{canopy}} \downarrow + \text{FAPAR}\_{\text{canopy}} \uparrow \tag{1}$$

where *FAPARcanopy*↓ describes the fraction of downwelling radiation absorbed within the canopy assuming the soil background is dark, and *FAPARcanopy*↑ describes the fraction of upwelling radiation absorbed within the canopy due to the interaction between the ground (soil and understory) and the canopy. A point worth emphasizing is that, assuming a black soil background, the *FAPARcanopy*↓ is equal to the *FAPARcanopy*.

Therefore, *FAPARcanopy*↓ can be given by:

$$FAPAR\_{\text{campy}} = (1 - \tau\_{PAI}) \ast (1 - Alledo\_{pure} \ast FVC) \tag{2}$$

where τ*PAI* is the transmittance of the whole canopy, which contains green and woody parts, *FVC* is the fraction of vegetation cover, and *Albedopure* is the visible (VIS) albedo for pure vegetation, which is represented by the VIS albedo for vegetation with a "saturated" LAI value (e.g., LAI ≥ 6) [14]. According to the statistical results by Liu et al. [15], *Albedopure* was set to 0.025 (for white-sky conditions) and 0.020 (for black-sky conditions) for all woody vegetation types (including evergreen needleleaf forest, evergreen broadleaf forest, deciduous needleleaf forest, and deciduous broadleaf forest).

τ*PAI* can be calculated as the product of τ*LAI* and τ*WAI*, and the directional transmittance can be determined using the gap fraction model [16,53,54]:

$$
\pi\_{\text{PAI}} = \pi\_{\text{LAI}} \times \pi\_{\text{WAI}} \tag{3}
$$

$$
\pi\_{LAI} = \text{\(e}^{-\text{kj}\times\text{G}(\text{\(0\)}\times\text{CI}\times\text{LAI}/\cos\left(\text{\(0\)}\right)}\tag{4}
$$

$$
\tau\_{WAI} = \varepsilon^{-k\_2 \times G(\theta) \times \mathbb{C}l \times \mathbb{W}Al / \cos \left( \theta \right)} \tag{5}
$$

where *LAI* and *WAI* are the leaf area index and woody area index, respectively, *k*1 and *k*2 are extinction coefficients for the green components and woody components, respectively, and θ is the solar zenith angle. *k*1 can be determined using the leaf absorptance in the VIS band and was set to 0.88 based on simulations made using the PROSPECT-5 model and also the Leaf optical properties experiment 93 (LOPEX'93) and Leaf Optical Properties Database (an experiment conducted at the National Institute for Agricultural Research in Angers, France in June 2003) [55–57]. The woody components were assumed to be opaque with a constant extinction coefficient (*k*2) of 0.91 based on the simulated LESS data. *G*(θ) is the projection of the unit foliage area on the plane perpendicular to the solar incident direction, θ. For green leaves, *G*(θ) is normally given a value of 0.5 for canopies with a spherical leaf angle distribution. For woody components, we assuming the woody components have the same angular distribution as green leaves, which also assumed by Chen et al. [58], Kucharik et al. [59], and Sea et al. [60]. CI is the clumping index; we also assume the same CI values for both green leaves and woody components based on the findings of Chen et al. [61] and Zou et al. [62].

The gap fraction model is also used to calculate FVC for a fixed solar zenith angle of 0◦ and a fixed value of 1 for the leaf extinction coefficient:

$$FVC \;=\; 1 - e^{-G(\theta) \times LAI \times CI} \tag{6}$$

Similarly, *FAPARcanopy*↑ can be calculated as follows:

$$FAPAR\_{\text{canopy}}\uparrow = (1 - \tau\_{PAI}) \ast (1 - Alledo\_{\text{pure}} \ast \text{FVC}) \ast \ \tau\_{PAI}^{\text{us}} \ast Alledo\_{\text{soil}} \tag{7}$$

where *Albedosoil* is the soil visible albedo, which can be generated by the NSM model proposed by Liu et al. [15]. τ*ws PAI* is the canopy transmittance under white-sky conditions, which can be calculated as

the product of the white-sky transmittance of leaves (τ*ws LAI*) and the white-sky transmittance of woody components (τ*ws WAI*):

$$\tau\_{LAI}^{\text{us}} = 2 \times \int\_0^{\frac{\pi}{2}} \left( e^{-k\_l \times G(\theta) \times LAI \times \Box / \cos(\theta)} \right) \times \sin(\theta) \times \cos(\theta) d\theta \tag{8}$$

$$\pi\_{\rm WAI}^{\rm ux} = 2 \times \int\_0^{\frac{\pi}{2}} \left( e^{-k\_2 \times G(\theta) \times \text{WAI} \times \text{CI} / \cos(\theta)} \right) \times \sin(\theta) \times \cos(\theta) d\theta \tag{9}$$

$$
\pi\_{PAI}^{\text{us}} = \pi\_{LAI}^{\text{us}} \times \pi\_{WAI}^{\text{us}} \tag{10}
$$

#### 2.3.2. Determination of Woody Area Index

In general, it is expensive and time-consuming to make accurate estimates of the woody area index (WAI), and destructive sampling is often the only option available for the quantification of the WAI in tropical evergreen forests [63]. Therefore, for generating global *FAPARgreen* datasets, the use of accurate estimates of WAI is unrealistic. Hence, in this study, we aimed to determine the global WAI using MCD15A2H LAI data and assumed that the WAI was constant within a given year.

First, we assumed that the woody-to-total area ratio is measured during the peak growth stage (July–August–September, JAS) when the LAI has its maximum value for the year. The woody-tototal area ratio was determined according to the forest types listed in the MODIS land cover product (MCD12Q1), including evergreen needleleaf forest (ENF), evergreen broadleaf forest (EBF), deciduous broadleaf forest (DBF), deciduous needleleaf forest (DNF) and mixed forest (MF). Therefore, WAI values for different forest vegetation types were calculated using a simple linear relationship between the plant area index (PAI) and the WAI:

$$PAI \;= \; PAI \times ratio\_{woody} \;= \; \frac{LAI\_{mw}}{1 - ratio\_{woody}} \times ratio\_{woody} \tag{11}$$

where *ratiowoody* is the mean value of the woody-to-total area ratio for various forest types, as given in the literature [58,64,65]. *LAImax* is the maximum value of LAI within a given year; this was acquired from MCD15A2H LAI products.

#### 2.3.3. Separating *FAPARgreen* and *FAPARwoody* from *FAPARcanopy*

In a similar way to Equation (1), the fraction of PAR absorbed by the green and woody components can also be separated into the fraction of PAR absorbed by the canopy for the downwelling radiation and that for the upwelling radiation reflected from the soil background. Firstly *FAPARcanopy*↓ is split into *FAPARgreen*↓ and *FAPARwoody*↓:

$$FAPAR\_{\text{canopy}\downarrow} = FAPAR\_{\text{gren}\downarrow} + FAPAR\_{\text{wody}\downarrow} \tag{12}$$

where *FAPARgreen*↓ and *FAPARwoody*↓ can be obtained as

$$FAPAR\_{\text{gren}\downarrow} = FAPAR\_{\text{campy}\downarrow} \times ratio\_{\text{gren}} \times w\_{1\downarrow} \tag{13}$$

$$FAPAR\_{\text{rowody}\downarrow} = FAPAR\_{\text{canopy}\downarrow} \times ratio\_{\text{rowody}} \times w\_{2\downarrow} \tag{14}$$

where *ratiogreen* and *ratiowoody* are the ratio of leaf area index to plant area index and woody area index to plant area index, respectively. *w*1↓ and *w*2↓ are the weighting coefficients for the green (i.e., photosynthetic) and woody components in terms of the radiation transfer within the canopy. The terms involved in equations (10) and (11) can be given as:

$$ratio\_{green} = \begin{array}{c} \frac{LAI}{PAI} \end{array} \tag{15}$$

$$ratio\_{\text{rowdy}} = \frac{\text{VAL}}{\text{PAI}} \tag{16}$$

$$
\pi\_{2\downarrow} = w\_{1\downarrow} \times \tau\_{LAI} \tag{17}
$$

$$w\_{1\downarrow} \times ratio\_{\text{green}} + w\_{2\downarrow} \times ratio\_{\text{woody}} = \mathbf{1} \tag{18}$$

*FAPARgreen*↓ and *FAPARwoody*↓ can then be obtained by solving equations (13)–(18):

$$FAPAR\_{\text{grav}\downarrow} = \frac{\text{ratio}\_{\text{grav}} \times FAPAR\_{\text{canopy}\downarrow}}{\text{ratio}\_{\text{grav}} + \text{τ}\_{LA} \times \text{τ}\_{\text{today}\_{\text{unody}}}} \tag{19}$$

$$FAPAR\_{\text{nvody}\downarrow} = \frac{\text{ratio}\_{\text{uvody}} \times FAPAR\_{\text{nvays}\downarrow} \times \tau\_{LAI}}{\text{ratio}\_{\text{gurn}} + \tau\_{LAI} \times \text{ratio}\_{\text{uvody}}} \tag{20}$$

Similarly, *FAPARgreen*↑ and *FAPARwoody*↑ can also be acquired:

$$FAPAR\_{\text{gren}\uparrow} = \frac{\text{ratio}\_{\text{gren}} \times \text{FAPAR}\_{\text{camppy}\uparrow} \times \text{trVAL}}{\text{ratio}\_{\text{wody}} + \text{rtOM} \times \text{ratio}\_{\text{gren}}} \tag{21}$$

$$FAPAR\_{\text{uvody}} = \frac{\text{ratio}\_{\text{uvody}} \times FAPAR\_{\text{annyp}}}{\text{ratio}\_{\text{uvody}} + \tau\_{\text{VA}\lambda} \times \text{ratio}\_{\text{gram}}} \tag{22}$$

Finally, *FAPARgreen* and *FAPARwoody* can be calculated as:

$$\text{FAPAR}\_{\text{green}} = \text{FAPAR}\_{\text{green}\uparrow} + \text{FAPAR}\_{\text{green}\downarrow} \tag{23}$$

$$FAPAR\_{\text{wody}} = FAPAR\_{\text{wody}} + FAPAR\_{\text{wody}} \tag{24}$$

#### **3. Results**

#### *3.1. Validation of the TriLay Method using Simulations made by the LESS Model*

It is too challenging to obtain in-situ measurements of *FAPARgreen* for forests, and so the simulated dataset (Table 1) derived using the LESS model was used for the validation of the TriLay model. Seventy-two black-sky simulations of *FAPARcanopy*, *FAPARgreen*, and *FAPARwoody* together with eight white-sky simulations were available for validation.

Figure 4 illustrates the validation results for *FAPARcanopy*, *FAPARgreen*, and *FAPARwoody*. The estimated and simulated FAPAR values are distributed close to the 1:1 line. Also, it can be seen that the TriLay approach can produce accurate estimates of *FAPARcanopy*, giving Root Mean Square Error (RMSEs) of 0.048 and 0.024 for black-sky and white-sky conditions, respectively, as against the LESS-simulated values. The corresponding R2 values are 0.945 and 0.999. For *FAPARgreen*, the validation results give RMSEs of 0.064 and 0.025, respectively, for black-sky and white-sky FAPAR. Finally, it can be seen that *FAPARwoody* can also be accurately estimated: the RMSE and R2 values are 0.042 and 0.709, respectively, for black-sky conditions, and 0.014 and 0.992 for white-sky conditions. These results show that the TriLay approach can be used to accurately estimate *FAPARcanopy*, *FAPARgreen*, and *FAPARwoody* for forest land cover types.

Furthermore, as illustrated in Figure 4b, there is a slight underestimation for *FAPARgreen* at smaller SZAs (0◦–60◦) and a slight overestimation for larger SZAs (70◦–80◦). From Figure 4c, it can be seen that *FAPARwoody* is also slightly underestimated at smaller SZAs (0◦–40◦) and overestimated at larger angles (50◦–80◦).

**Figure 4.** Validation of *FAPARcanopy*, *FAPARgreen* and *FAPARwoody* estimates made by the Trilay model against the LargE-Scale remote sensing data and image simulation (LESS)—simulated FAPARs: (a–c) black-sky *FAPARcanopy*, *FAPARgreen*, and *FAPARwoody*; (d–f) white-sky *FAPARcanopy*, *FAPARgreen*, and *FAPARwoody*.

#### *3.2. Comparison of Di*ff*erent Methods using the LESS Simulations*

Hall et al. [38] estimated *FAPARgreen* based on a simple linear relationship between *FAPARcanopy* and *FAPARgreen*. They used the ratio *LAIgreen*/*LAItotal* to determine *FAPARgreen*:

$$FAPAR\_{green} = FAPAR\_{canopy} \times \frac{LAI\_{green}}{total\ LA} \tag{25}$$

where total LAI is the PAI mentioned above. In order to test the accuracy of the linear mixture method, we also validated the *FAPARgreen* estimated by the linear mixture method [38] using LESS-simulated data.

Figure 5 shows the accuracy assessment results for the linear mixture method. The results show a noticeable underestimation for *FAPARgreen* and a huge overestimation for *FAPARwoody*. Although the FAPARs retrieved using the linear mixture method are highly correlated with the LESS simulation, with R2 values of 0.979, 0.996 for *FAPARgreen*, and 0.934, 0.985 for *FAPARwoody* under black-sky and white-sky conditions, the corresponding RMSE values are much higher than those found using our TriLay approach—0.114, 0.106 for *FAPARgreen* as against 0.064, 0.025, and 0.113, 0.106 for *FAPARwoody* as against 0.042, 0.014 under the black-sky and white-sky conditions, respectively.

**Figure 5.** Validation of *FAPARgreen* and *FAPARwoody* estimated by the linear method using the LESS-simulated values of FAPAR: (a–b) black-sky *FAPARgreen* and *FAPARwoody* against the LESS-simulated *FAPARgreen*; (c–d) white-sky *FAPARgreen* and *FAPARwoody* against the LESS-simulated *FAPARwoody*.

The FAPAR values derived without considering the woody components (*FAPARnoWAI*) were also validated using the LESS simulations; the results are shown in Figure 6. These results give an RMSE of 0.071 and R2 of 0.920 for the black-sky conditions and the corresponding values of 0.043 and 0.999 for the white-sky conditions. Although the R2 values are higher, the RMSEs are still greater than those found using the TriLay approach.

**Figure 6.** Validation of *FAPARnoWAI* using the LESS-simulated FAPAR values: (**a**) black-sky and (**b**) white-sky *FAPARnoWAI* against the LESS-simulated *FAPARgreen*.

In order to further compare the performances of different methods, we also calculated the bias for *FAPARgreen*, *FAPARwoody*, and *FAPARnoWAI* using the LESS-simulated FAPAR values. Table 2 summarizes the retrieval accuracy for these three methods. The results show that the TriLay method gave the best FAPAR retrieval results, having the smallest bias and RMSE values (RMSE=0.064 and 0.025, and bias = −6.02%, −4.04% for *FAPARgreen* under black-sky and white-sky conditions, respectively).

**Table 2.** Retrieval accuracy of *FAPARgreen*, *FAPARwoody*, and *FAPARnoWAI* validated using the LESS simulations.

#### *3.3. Temporal Variations in Di*ff*erent FAPAR Products*

Using the TriLay approach, black-sky and white-sky products, including *FAPARgreen*, *FAPARwoody*, and *FAPARcanopy*, and also FAPAR without woody components (*FAPARnoWAI*), were generated for tile h10v05 in 2017. The black-sky FAPAR products were determined according to the SZA values at 10:30 am local time (the overpass time of the Terra satellite).

To investigate the seasonal variations in *FAPARgreen* and the other FAPAR products, only the mean black-sky FAPAR values were calculated for the different forest vegetation types, as shown in Figure 7. The proportion of black-sky *FAPARgreen* in *FAPARcanopy* and *FAPARnoWAI*, and also the bias between the black-sky *FAPARgreen*, *FAPARcanopy*, and *FAPARnoWAI* were calculated for different periods during 2017 in order to analyze the differences between different black-sky FAPAR products, as well as to quantify the contribution of the woody components for several forest vegetation types (deciduous broadleaf forest, deciduous needleleaf forest, evergreen broadleaf forest, and evergreen needleleaf forest, referred to as DBF, DNF, EBF, and ENF, respectively). The results are illustrated in Figure 8 and Table 3.

In general, the black-sky *FAPARgreen* and *FAPARnoWAI* exhibit typical seasonal variations for the selected forest types, with low values during the early and late period (January–February–March (JFM) and October–November–December (OND)) and high values during the peak growth stage (July–August–September (JAS)). The black-sky *FAPARcanopy* is obviously higher than the black-sky *FAPARgreen* during the whole year and has a much smaller seasonal variation. The black-sky *FAPARwoody* behaves in the opposite way; for deciduous seasons (JFM and OND), *FAPARwoody* is close to or higher than its value during the peak growth stage (JAS) because the black-sky FAPAR increases with increasing SZA value. (The mean SZA at 10:30 am varies from 22.98◦ (22 December) to 62.25◦ (22 June) within tile h10v05.)

The black-sky *FAPARgreen* is about 52.59% and 60.60% of the black-sky *FAPARcanopy* for deciduous broadleaf forest during JFM and OND, respectively; the corresponding figures for deciduous needleleaf forest are only about 23.86% and 35.75%. During the peak growth stage (JAS), the black-sky *FAPARgreen* is about 93.36%, 75.02%, 90.93% and 87.14% of *FAPARcanopy* for DBF, DNF, EBF, and ENF, respectively.

There are also small discrepancies between *FAPARnoWAI* and *FAPARgreen* (Figure 8). In particular, for deciduous needleleaf forests, the black-sky *FAPARnoWAI* is overestimated by 38.30% and 28.46% during the early and late stages of the year (JFM and OND). For evergreen forests, the difference can be neglected as there is only a very slight underestimation of 0.68% to 2.39% during the whole year.

**Figure 7.** Temporal variations in the mean black-sky *FAPARgreen*, *FAPARnoWAI*, *FAPARcanopy*, and *FAPARwoody* products for forest vegetation types within tile h10v05 during 2017 (tile h10v05 is located in North America, covering 30.00◦ N–40.00◦ N and 80.00◦ W–104.4◦ W).

**Figure 8.** Bias between the black-sky *FAPARgreen* and other black-sky FAPAR products within tile h10v05 during different periods in 2017: (**a**) bias between *FAPARgreen* and *FAPARnoWAI*; (**b**) bias between *FAPARgreen* and *FAPARcanopy*. JFM, AMJ, JAS, and OND represent the four seasons January to March, April to June, July to September and October to December, respectively.

**Table 3.** The ratios of black-sky *FAPARgreen* to black-sky *FAPARcanopy* (*Rcanopy*) and to *FAPARnoWAI* (*RnoWAI*) for selected forest types for different periods of 2017.

JFM, AMJ, JAS, and OND represent the four seasons January to March, April to June, July to September and October to December, respectively.

#### **4. Discussion**

#### *4.1. Uncertainty in Determining WAI*

Currently, methods for measuring the woody area index (WAI) and woody-to-total area ratio can be classified into direct methods (e.g., destructive sampling) and indirect methods [66]. However, both the direct and indirect methods can only be applied at small scales and to a limited range of vegetation types. In this study, we used a constant value of the woody-to-total area ratio for each forest type and derived the WAI from the PAI value for the peak growth stage. We obtained woody-to-total area ratios for different forest types (ENF, EBF, DBF, and DNF) from an extensive literature review [58,64–67]—the statistical metrics, including the mean, standard deviation, and coefficient of variation, are shown in Table 4. These results show that the woody-to-total area ratios for different forest types vary from 0.158 to 0.3. The variation in this ratio within each of the forest types is also large—10.00% to 82.22%. Therefore, the woody-to-total area ratio not only varies with the forest type but also changes a lot for each individual forest type. This means that there is definitely some uncertainty due to these factors. Even so, Zou et al. [66] showed that the woody-to-total area ratio is relatively stable for the same forest stand and thus the assumption of a fixed woody-to-total area ratio for each forest type is reasonable. Differences in in situ measurement methods can also contribute to the variation in the ratio for a given forest type. The retrieval of *FAPARgreen* can thus be improved if an accurate woody-to-total area ratio dataset is available.

**Table 4.** Statistical details of prior woody-to-total area ratios.

#### *4.2. Uncertainty Caused by the Use of Fixed Values of the Extinction Coe*ffi*cients and Albedopure*

The canopy directional transmittance can be determined using the gap fraction model [16,53,54]. Splitting the canopy into green components and woody components requires that the transmittances are also calculated separately. For the green components, the leaf absorptance at the VIS band varies from 0.79 to 0.94 as the chlorophyll content varies (from 20 to 100 μg/cm2 under natural conditions), according to figures obtained using the PROSPECT-5 model [56,57]. Also, because no remote sensing leaf chlorophyll content product is available, the extinction coefficient for leaves (*k*1) was assumed to have a fixed value of 0.88 (corresponding to a chlorophyll content of 35 μg/cm2 [15]). For the extinction coefficient of woody components (*k*2), we also used a fixed value of 0.91, based on the data simulated by LESS (a needleleaf forest scene). However, Suwa et al. [68] reported a smaller extinction coefficient of 0.77 for brighter woody stems. Therefore, the extinction coefficient of woody components may vary with forest type, and thus the use of a fixed value for *k*2 is also a source of error. At present, it is still very challenging to determine extinction coefficients for different forest canopy types.

In addition, a fixed value of *Albedopure* was used in the TriLay model for woody vegetation types (i.e., ENF, EBF, DBF, and DNF), and was approximated based on the dense vegetation [15]. According to the statistical results obtained from the MCD43A3 albedo product by Liu et al. [15], the visible albedo of dense vegetation with a "saturated" LAI value (e.g., LAI = 6) is very low and stable, with a mean value of 0.025 and a variance of 0.007 for white-sky condition, and a mean value of 0.020 and a variance of 0.006 for black-sky conditions. Therefore, the use of the prior VIS albedo values for "pure" vegetation may introduce a very small error, but it should be negligible [14].

#### *4.3. Setting the Clumping Index for Photosynthetic and Woody Components*

The clumping index characterizes the grouping of foliage within distinct canopy structures (such as tree crowns, shrubs, and row crops) relative to a random spatial distribution of leaves and is an important structural parameter for plant canopies that can influence canopy radiation regimes [49]. In our TriLay approach, the clumping index for woody components was assumed to be the same as for green leaves, and so is another definite source of error [69]. However, it is currently challenging to obtain the clumping effects of woody components within forest canopies [62]. Furthermore, Chen et al. [61] indicated that the clumping of shoots in branches has a similar effect to the clumping of leaves within shoots. What's more, Zou et al. [62] found that the differences between the estimated CI for canopy and woody components was below 6% at the zenithal ranges of 0◦–75◦, and the difference was only 2% in the range of 30◦–60◦, which is quite small at most medium zenithal ranges thus can represent most actual conditions. Based on these results and the unavailability of CI datasets for forest canopies at large scales, we directly used the value of CI from Jiao et al. [49] to describe the clumping effect for both leaves and woody components. Therefore, the use of the same clumping index for both leaves and woody components is reasonable and credibe enough.

#### **5. Conclusions**

In this paper, a triple-source leaf–wood–soil layer (TriLay) method for separating *FAPARgreen* and *FAPARwoody* from *FAPARcanopy* using the MODIS LAI, land cover, and non-linear spectral mixture model (NSM)-retrieved soil albedo [15] together with global CI products [49] was proposed.

According to the validation carried out using LESS-simulated FAPAR values, the TriLay *FAPARgreen* was more accurate (R2 = 0.937, RMSE = 0.064 and bias = −6.02% for black-sky conditions; R2 = 0.997, RMSE = 0.025 and bias = −4.04% for white-sky conditions) than the traditional linear method (R2 = 0.979, RMSE = 0.114 and bias = −18.04% for black-sky conditions; R2 = 0.996, RMSE = 0.106 and bias = −16.93% for white-sky conditions), and also more accurate than FAPAR obtained without the consideration of woody components (*FAPARnoWAI*) (R2 = 0.920, RMSE = 0.071 and bias = −7.14% for black-sky conditions; R2 = 0.999, RMSE = 0.043 and bias = −6.41% for white-sky conditions). A comparison of the results for black-sky *FAPARgreen* against *FAPARnoWAI* and *FAPARcanopy* showed that the discrepancies between the black-sky *FAPARgreen* and other FAPAR products could not be ignored for forest types. In particular, for deciduous needleleaf forest, the black-sky *FAPARgreen* contributed only about 23.86% and 35.75% of *FAPARcanopy* during the early and late stages (JFM and OND) of the year, respectively, and 75.02% during the peak growth stage (JAS). There were also smaller discrepancies between the black-sky *FAPARnoWAI* and *FAPARgreen*. For deciduous needleleaf forests, in particular, the black-sky *FAPARnoWAI* was overestimated by 38.30% and 28.46%, respectively, during the early and late stages of the year (JFM and OND).

Overall, this study provides a new method for partitioning *FAPARcanopy* into *FAPARgreen* and *FAPARwoody* for forest types and will improve the understanding of energy exchange within the canopy. In addition, the exclusion of the contribution of woody components may certainly improve the accuracy of the *FAPARgreen* estimates for forest types, which is significant in terms of the better modeling of vegetation photosynthesis.

**Author Contributions:** Conceptualization—S.C. and L.L.; methodology—S.C. and L.L.; software—S.C., X.Z. and X.C.; validation—S.C., Y.X. and D.X.; formal analysis—S.C. and X.Q.; investigation—S.C. and L.L.; resources—S.C.; writing—original draft preparation—S.C.; writing—review and editing—L.L. and X.L.

**Funding:** This research was funded by the National Key Research and Development Program of China, grant number 2017YFA0603001, the Key Research Program of the Chinese Academy of Sciences (ZDRW-ZS-2019-1), the National Natural Science Foundation of China (41825002).

**Acknowledgments:** The authors gratefully acknowledge the CI products provided by Ziti Jiao at Beijing Normal University, and the MODIS products obtained from the Land Products Land Processes Distributed Active Archive Center (LP DAAC).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Climate Data Records of Vegetation Variables from Geostationary SEVIRI**/**MSG Data: Products, Algorithms and Applications**

**Francisco Javier García-Haro 1,\*, Fernando Camacho 2, Beatriz Martínez 1, Manuel Campos-Taberner 1, Beatriz Fuster 2, Jorge Sánchez-Zapero 2 and María Amparo Gilabert 1**

Received: 26 July 2019; Accepted: 6 September 2019; Published: 9 September 2019

**Abstract:** The scientific community requires long-term data records with well-characterized uncertainty and suitable for modeling terrestrial ecosystems and energy cycles at regional and global scales. This paper presents the methodology currently developed in EUMETSAT within its Satellite Application Facility for Land Surface Analysis (LSA SAF) to generate biophysical variables from the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) on board MSG 1-4 (Meteosat 8-11) geostationary satellites. Using this methodology, the LSA SAF generates and disseminates at a time a suite of vegetation products, such as the leaf area index (LAI), the fraction of the photosynthetically active radiation absorbed by vegetation (FAPAR) and the fractional vegetation cover (FVC), for the whole Meteosat disk at two temporal frequencies, daily and 10-days. The FVC algorithm relies on a novel stochastic spectral mixture model which addresses the variability of soils and vegetation types using statistical distributions whereas the LAI and FAPAR algorithms use statistical relationships general enough for global applications. An overview of the LSA SAF SEVIRI/MSG vegetation products, including expert knowledge and quality assessment of its internal consistency is provided. The climate data record (CDR) is freely available in the LSA SAF, offering more than fifteen years (2004-present) of homogeneous time series required for climate and environmental applications. The high frequency and good temporal continuity of SEVIRI products addresses the needs of near-real-time users and are also suitable for long-term monitoring of land surface variables. The study also evaluates the potential of the SEVIRI/MSG vegetation products for environmental applications, spanning from accurate monitoring of vegetation cycles to resolving long-term changes of vegetation.

**Keywords:** meteosat second generation (MSG); biophysical parameters (LAI; FVC; FAPAR); SEVIRI; climate data records (CDR); stochastic spectral mixture model (SSMM); Satellite Application Facility for Land Surface Analysis (LSA SAF)

#### **1. Introduction**

Observations and monitoring are essential components of Global Framework for Climate Services initiative [1,2] that contribute significantly to meet the needs of climate research services. The role that satellites have played in observing the variability and change of the Earth system has increased significantly in the last decades. Thanks to the Earth observation (EO) satellite systems, a large number of variables related to the atmosphere, oceanic and terrestrial domains are accessible [3].

The Satellite Application Facility (SAF) for Land Surface Analysis (LSA) is a leading center for development and operational retrieval of land surface variables from European Organization for the Exploitation of Meteorological Satellites (EUMETSAT), specifically aimed to understand and quantify terrestrial processes and land-atmosphere interactions. Land surface models' assessment and improvement of high quality and robust EO products are among the target applications of LSA SAF products. Particularly, the geostationary Spinning Enhanced Visible and Infrared Imager (SEVIRI) on board the Meteosat Second Generation (MSG) platform provides a very high temporal resolution (15 minutes) required for resolving the diurnal cycles of radiation, surface albedo and surface temperature [4]. The MSG program covers a series of four identical satellites MSG 1-4 (Meteosat 8-11) launched in 2002, 2005, 2012, and 2015 respectively. SEVIRI operates in 12 spectral channels [5], although only channels 1 (VIS0.6), 2 (VIS0.8) and 3 (NIR1.6), centered at about 0.635 μm (red), 0.81 μm (NIR) and 1.64 μm (SWIR), respectively situate in the optical domain.

The LSA SAF provides also variables for the characterization of terrestrial ecosystems and their role in the energy balance of Earth, such as land surface fluxes and three important biophysical variables related with the amount, structure and state of vegetation: Leaf Area Index (LAI), Fraction of Absorbed Photosynthetically Active Radiation (FAPAR) and Fractional Vegetation Cover (FVC) [6]. These variables are widely used in many land-surface vegetation, climate, and crop production models [7–13]. FVC represents the fraction of green vegetation covering a unit area of horizontal soil, corresponding to the gap fraction in the nadir direction. LAI is a measure of the amount of live foliage present in the canopy per unit ground surface. The FAPAR is the fraction of photosynthetically active radiation (400-700 nm) absorbed by the green parts of the canopy, and therefore constitutes an indicator of the presence, health and productivity of live vegetation.

Nowadays, consistent long-term data records of these variables with well-characterized uncertainty are currently required by the scientific community to model terrestrial ecosystems and energy cycles at regional and global scales [14]. Models related with land surface and climate usually require Near-Real-Time (NRT) datasets of these variables at high temporal frequencies, and with regional to global coverage [15]. Both, LAI and FAPAR are recognized as essential climate variables in the terrestrial domain by Global Climate Observing System (GCOS) [16].

Several methods have been proposed for estimating biophysical parameters from large scale optical sensors, some of which have been implemented in operational processing lines from POLarization and Directionality of Earth Reflectances (POLDER) [17], Moderate Resolution Imaging Spectroradiometer (MODIS) [18], MISR (Multi-angle Imaging SpectroRadiometer) [19], MEdium Resolution Imaging Spectrometer (MERIS) [20,21], Sea-viewing Wide Field-of-view *Sensor* (SEAWIFS) [22], VEGETATION [23,24] and Advanced Very High Resolution Radiometer (AVHRR) [25]. These methods can be roughly classified in two main categories: (i) methods based on optimized vegetation indices [26], and (ii) methods based on the inversion of physical models. Both approaches present advantages and drawbacks. It is referred to [27] for a detail intercomparison on the spread of currently available physical models. Although methods that invert physical models can potentially be applied to varying surface types, they require specifying a large number of model parameters, which may lead to unstable solutions [28]. Most inversion methods use one-dimensional (turbid) models (e.g., PROSAIL [29]), thereby assuming surface homogeneity within an image pixel. They are thus limited to address the sub-pixel mixing of land cover types, which is common at coarse resolutions satellites. The process of identifying the sub-pixel proportions of the constituent components is called spectral mixture model (SMM). SMM approaches are methods especially adequate for global studies, since the spatial variability within pixel is high [30,31].

In this paper, we describe the methodology currently used in LSA SAF to operationally generate the suite of SEVIRI/MSG (FVC, LAI and FAPAR) vegetation products. The high frequency of acquisition provided by the SEVIRI instrument assures the availability of cloud-free data for appropriately monitoring the vegetation dynamics in large regions. When one is choosing a retrieval method, one needs to consider the limited spatial and spectral resolution of SEVIRI instrument (i.e., only 3 optical channels and a spatial resolution of 3 km at the sub-satellite point that declines rapidly at higher latitudes). Furthermore, the challenges that the MSG geostationary satellite poses as compared to polar-orbiting satellite sensors due to diurnal variations in the solar zenith angle need to be considered. The algorithm to retrieve FVC relies on a novel stochastic SMM method, which addresses the large variability of soils and vegetation types using statistical (sum of Gaussian) distributions. The LSA SAF

operational scheme to retrieve LAI and FAPAR uses statistical relationships general enough for global applications, proposed by Roujean and Lacaze [16] and Roujean and Bréon [32], respectively.

LSA SAF has been providing in NRT vegetation products (FVC, LAI and FAPAR) at two different time resolutions (daily and 10-days) since 2008 over the geostationary Meteosat disk, covering Europe, Africa, the Middle East and parts of South America. The suite of SEVIRI/MSG vegetation products is freely disseminated to users through LSA SAF website (https://landsaf.ipma.pt). During the different operational phases, these vegetation products have been routinely assessed and some improvements have been made in the upgraded versions. LSA SAF has recently reprocessed the entire MSG archive with the latest version of the several retrieval algorithms in order to obtain a continuous and homogeneous Climate Data Records (CDR) of vegetation products. More than fifteen years (January 2004- present) of homogeneous consistent estimates of FVC, LAI and FAPAR fields are available. The 10-days (FVC, LAI and FAPAR) CDRs are part of the suite of EUMETSAT climate products, which include Fundamental CDRs, Thematic CDRs and Data for operational climate monitoring and are also available from EUMETSAT web page (https://navigator.eumetsat.int/).

The paper presents the following structure. Section 2 is dedicated to describe the algorithms to retrieve FVC, LAI and FAPAR. The suite of LSA SAF SEVIRI/MSG vegetation products are described in Section 3, including expert knowledge and a quality assessment of the consistency between the products. Section 4 provides further insight about the added value and introduces examples of its application in numerous environmental applications.

#### **2. Algorithm Description**

ε

#### *2.1. SEVIRI*/*MSG*

The algorithm for retrieving biophysical variables from SEVIRI/MSG data uses as input normalized spectral bidirectional reflectance factor (BRDF) at a fixed angular configuration (i.e., vertical illumination and observation angles of 0◦) in order to reduce angular effects. This hotspot geometry leads to a minimum contribution of the shadow proportion and a precise estimation of FVC, coinciding with the complement to unity of the gap fraction at nadir direction [16]. Estimating the FVC with increased sun zenith angles would lead to an overestimation of FVC. The main drawback of this geometry is the significant contribution of illuminated soil background, which constitutes a confounding factor for the retrieval of FVC. Figure 1 outlines the main ingredients for deriving the SEVIRI/MSG products including the stochastic SMM approach for retrieving FVC, a pragmatic method for LAI retrieval from FVC, and a method based on an optimized vegetation index for FAPAR determination.

**Figure 1.** Flow chart of the algorithm for FVC, LAI and FAPAR determination.

ε

The BRDF algorithm applies a semi-empirical reflectance model in order to invert SEVIRI/MSG top-of-canopy reflectance factor values into three parameters (*k0*, *k1*, *k2*). The input of the FVC algorithm is atmospherically corrected BRDF *k0* parameters in the MSG channels 1, 2 and 3, while FAPAR uses all three BRDF parameters (*k0*, *k1*, *k2*) in order to correct surface's reflectance anisotropy and minimize the effect of soil reflectance (see Figure 1). The algorithm to correct atmospheric effects and generate BRDF model parameters is fully described in [33], and in the ATBD of the Albedo MSG product [34].

#### *2.2. FVC Algorithm*