Identiﬁcation of Robust Hybrid Inversion Models on the Crop Fraction of Absorbed Photosynthetically Active Radiation Using PROSAIL Model Simulated and Field Multispectral Data

: The fraction of absorbed photosynthetically active radiation (FPAR), which represents the capability of vegetation-absorbed solar radiation to accumulate organic matter, is a crucial indicator of photosynthesis and vegetation growth status. Although a simpliﬁed semi-empirical FPAR estimation model was easily obtained using vegetation indices (VIs), the sensitivity and robustness of VIs and the optimal inversion method need to be further evaluated and developed for canola FPAR retrieval. The objective of this study was to identify the robust hybrid inversion model for estimating the winter canola FPAR. A ﬁeld experiment with different sow dates and densities was conducted over two growing seasons to obtain canola FPARs. Moreover, 29 VIs, two machine learning algorithms and the PROSAIL model were incorporated to establish the FPAR inversion model. The results indicate that the OSAVI, WDRVI and mSR had better capability for revealing the variations of the FPAR. Three parameters of leaf area index (LAI), solar zenith angle (SZA) and average leaf inclination angle (ALA) accounted for over 95% of the total variance in the FPARs and OSAVI exhibited a greater resistance to changes in the leaf and canopy parameters of interest. The hybrid inversion model with an artiﬁcial neural network (ANN-VIs) performed the best for both datasets. The optimal hybrid inversion model of ANN-OSAVI achieved the highest performance for canola FPAR retrieval, with R 2 and RMSE values of 0.65 and 0.051, respectively. Finally, the work highlights the usefulness of the radiation transfer model (RTM) in quantifying the crop canopy FPAR and demonstrates the potential of hybrid model methods for retrieving the canola FPAR at each growth stage.


Introduction
The FPAR is the ratio of photosynthetically active radiation absorbed by the green part of the vegetation canopy to total photosynthetically active radiation (PAR) within the wavelength range of 400-700 nm. It is an important parameter for describing vegetation photosynthesis and physiological characteristics [1]. Since the 1990s, the FPAR has been extensively used in various areas, including assessing vegetation productivity [2], guiding vegetation growth management [3], and diagnosing drought and climate change [4]. Moreover, the FPAR has been recognized as one of the key climate variables that influence global climate change according to the Global Climate Observation System (GCOS) [5]. Therefore, accurate and reliable FPAR acquisition is of great significance for depicting terrestrial ecosystem processes and estimating crop yields [6].
FPAR estimation methods are usually divided into two categories: empirical methods based on vegetation indices (VIs) and physical model inversion based on the radiation transfer model (RTM). Vegetation indices can quantify the absorption and reflection of solar radiation by a vegetation canopy and are widely used for estimating the FPAR due to their fewer parameters and high computational efficiency [7]. Over the past few decades, numerous Vis, such as the normalized difference VI (NDVI) [8], enhanced VI (EVI) [9] and simple ratio index (SR) [10], have been developed, evaluated and shown to be effective in estimating the FPARs for various vegetation types. However, there are still several issues associated with the VI approach that need to be addressed. First, while VIs can reveal the physiological and ecological characteristics of vegetation, they cannot fully capture the changes in solar effective radiation within the canopy and surface [11,12]. Furthermore, estimating VIs introduced uncertainty due to factors like variations in solar elevation angle and soil background. Moreover, most VI methods were statistically based and exhibited different relationships between VIs and FPARs for various spatial and temporal vegetation dynamics, resulting in regional restrictions in their application. Glenn et al. [13] suggested that VIs could be simply used to determine canopy light absorption without involving the canopy structural level.
In contrast, the RTM could elucidate and simulate the process of FPAR change by incorporating explicit geometric models of the canopy structure [14]. Meng et al. [15] quantitatively retrieved the spatial distribution of the absorbed FPAR in salt marsh vegetation from satellite images using the soil-canopy-photochemical and energy flux observation coupling model (SCOPE). Liu et al. [16] employed the modified simultaneous heat and water radiation model (MSHAW) to simulate the vertical profile of FPAR within the canopy. The RTM demonstrated greater adaptability to diverse environments and background conditions for FPAR simulation, but the model uncertainties associated with the input parameters also affected the RTM accuracy. Computational technology advances, more optimization algorithms (e.g., gradient methods, genetic algorithms and simulated annealing) and machine learning algorithms (artificial neural networks (ANNs), support vector regression (SVR) and random forest (RF)) have been combined with RTMs to improve the FPAR inversion accuracy [17,18]. For example, Verger et al. [19] used the PROSAIL model in conjunction with a neural network algorithm to investigate a series of measured crop FPARs and model the simulated FPARs.
In addition to the mechanism-based inversion model, other studies also evaluated robust estimation models through sensitivity analysis, which investigates the impacts of parameter variations on simulation results (e.g., band reflectance and VI) [20,21]. For example, Leolini et al. [22] constructed the relationship between the rescaled VIs (NDVI, OSAVI, etc.) and the observed FPAR, and validated this relationship using independent datasets. Dong et al. [23] analyzed the sensitivity and correlation between the three modified VIs and the FPAR using reflectance data simulated by the PROSAIL model. However, most sensitivity analyses of VIs primarily focus on their influence on LAI changes, and limited studies on FPAR changes have been reported. Although LAI and FPAR are strongly correlated based on the Beer-Lambert law, the relationship is also modulated by the extinction coefficient, introducing a certain degree of uncertainty. Moreover, most previous studies focused on a single empirical or physical method, but few combined multiple methods for FPAR inversion. Hou et al. [24] used the PROSAIL model and a look-up table (LUT) algorithm to investigate the FPAR inversion from satellite data for different vegetation types in various regions. Kolassa et al. [25] proposed an empirical model calibration method to distinguish vegetation types based on the deviation between the simulated FPAR and MODIS-observed FPAR. Furthermore, studies on crop FPAR retrieval mainly focused on crops such as wheat, maize and soybean [26,27], but limited studies were performed on canola. Monitoring canola FPAR dynamics could effectively diagnose growth status and stage nodes, providing key variables for crop productivity models. Due to differences in the canopy structure between different crops, existing FPAR estimation models developed for other crops need further evaluation and validation when applied to canola [28,29].
Therefore, to identify the robust hybrid estimation methods for improving the canola FPAR retrieval accuracy, the PROSAIL model, which is a classical model mechanism of radiative transfer for vegetation canopy, was employed and recompiled to simulate FPAR variations. A total of 29 VIs and two machine learning algorithms were devoted to retrieving the FPAR. The model-estimated and in situ measured FPAR from a two-year variations. A total of 29 VIs and two machine learning algorithms were devoted to retrie ing the FPAR. The model-estimated and in situ measured FPAR from a two-year fie experiment were used to evaluate the inversion performance. The main objectives of t study were to (1) reveal the capability of different VIs in interpreting FPAR variation, analyze the sensitivity of VIs to changes in leaf and canopy parameters, and (3) ident the robust hybrid estimation methods that consider inversion accuracy in different can growth stages.

Study Area
From 2020-2022, a field experiment on canola was conducted at the Yangzijin E logical Experimental Field of Yangzhou University. The experimental station was locat in the middle of Jianghuai Plain, at coordinates 119°24′ E and 32°21′ N, with an elevati of 5 m above sea level and belonging to the subtropical monsoon climate zone (Figure The annual average frost-free period, air temperature, evaporation and precipitation we 223 days, 14.8 °C, 937 mm and 1063 mm, respectively. The soil type in the experimen field was sandy loam, with pH value, field capacity, wilting point and bulk density of 8. 25%, 7.8% and 1.24 g cm −3 , respectively. The initial nutrient contents in the top 20 cm soil surface were as follows: organic matter 10.2 g kg −1 , total nitrogen 0.97 g kg −1 , availa phosphorus 16.3 mg kg −1 and available potassium 151.2 mg kg −1 [30].

Experimental Design
In this experiment, the canola variety tested was "Deaiyou 558". The 12 different so ing scenarios were designed and performed in the winter canola growing season (Figu 2). Referring to local practices, four sowing dates and three sowing densities were creat to modulate the crop growth process. For the sowing date treatments, four periods d noted early (ES 21 September), normal (NS1 6 October, NS2 23 October) and late (LS November) were chosen. Regarding the sowing density treatments, low (LD 12.5 pla m −2 ), medium (MD 25 plants m −2 ) and high densities (HD 37.5 plants m −2 ) were adopte Thus, a total of 12 treatments and 36 plots (each treatment with three replicates) we established. Each plot area was 5 m × 5 m and the interval between adjacent plots wa m. The seeds were manually sown with a row spacing of 40 cm and a plant spacing of cm. After the emergence of canola seedlings, the seedlings were thinned at the 3-leaf sta

Experimental Design
In this experiment, the canola variety tested was "Deaiyou 558". The 12 different sowing scenarios were designed and performed in the winter canola growing season ( Figure 2). Referring to local practices, four sowing dates and three sowing densities were created to modulate the crop growth process. For the sowing date treatments, four periods denoted early (ES 21 September), normal (NS1 6 October, NS2 23 October) and late (LS 6 November) were chosen. Regarding the sowing density treatments, low (LD 12.5 plants m −2 ), medium (MD 25 plants m −2 ) and high densities (HD 37.5 plants m −2 ) were adopted. Thus, a total of 12 treatments and 36 plots (each treatment with three replicates) were established. Each plot area was 5 m × 5 m and the interval between adjacent plots was 2 m. The seeds were manually sown with a row spacing of 40 cm and a plant spacing of 10 cm. After the emergence of canola seedlings, the seedlings were Agronomy 2023, 13, 2147 4 of 22 thinned at the 3-leaf stage and the seedlings were fixed at the 5-leaf stage. Other field management practices, such as weeding, insecticide and disease prevention, followed the local field practices. nomy 2023, 13, x FOR PEER REVIEW 4 of and the seedlings were fixed at the 5-leaf stage. Other field management practices, su as weeding, insecticide and disease prevention, followed the local field practices.

Canopy Multispectral Image Acquisition
During the growth period of the canola, multispectral images were obtained usi the Rededge-mx multi-spectrometer (MicaSense, Inc., Seattle, WA, USA), which w mounted on a DJI Inspire 2 UAV remote sensing platform (DJI Inc., Shenzhen, Chin ( Figure 3). The spectral images were captured approximately 7-10 days after sowing, an the UAV remote sensing operation was carried out from 11 am to 1 pm, specifically sunny and windless weather conditions. During each flight, five images with a resoluti of 1280 × 960 pixels were obtained, corresponding to the blue (475 nm), green (560 nm red (668 nm), red edge (717 nm) and near-infrared (840 nm) bands. During each flig operation, images of a reference plate were synchronously captured for radiation calibr tion. The flight height during the operation was set at 20 m, with an along-and cross-tra overlap of 75%. The camera lens was pointed in the downward direction, capturing mu tispectral images with a ground resolution of 1.4 cm at a field of view of 47.2°. After ea flight, the raw images were automatically processed using Pix4Dmapper software, inclu ing camera optimizing, photo aligning, radiometric correction and orthomosaic buildin Finally, image reflectance data for each pixel throughout crop growth seasons were ge erated. For each plot, the average band reflectance was extracted using the selected regi of interest (ROI) that excluded two rows of canola from the boundary area to avoid inte ference from plot edge effects, where crops grew more vigorously due to better lightin ventilation and population competition.

Canopy Multispectral Image Acquisition
During the growth period of the canola, multispectral images were obtained using the Rededge-mx multi-spectrometer (MicaSense, Inc., Seattle, WA, USA), which was mounted on a DJI Inspire 2 UAV remote sensing platform (DJI Inc., Shenzhen, China) ( Figure 3). The spectral images were captured approximately 7-10 days after sowing, and the UAV remote sensing operation was carried out from 11 am to 1 pm, specifically in sunny and windless weather conditions. During each flight, five images with a resolution of 1280 × 960 pixels were obtained, corresponding to the blue (475 nm), green (560 nm), red (668 nm), red edge (717 nm) and near-infrared (840 nm) bands. During each flight operation, images of a reference plate were synchronously captured for radiation calibration. The flight height during the operation was set at 20 m, with an along-and cross-track overlap of 75%. The camera lens was pointed in the downward direction, capturing multispectral images with a ground resolution of 1.4 cm at a field of view of 47.2 • . After each flight, the raw images were automatically processed using Pix4Dmapper software, including camera optimizing, photo aligning, radiometric correction and orthomosaic building. Finally, image reflectance data for each pixel throughout crop growth seasons were generated. For each plot, the average band reflectance was extracted using the selected region of interest (ROI) that excluded two rows of canola from the boundary area to avoid interference from plot edge effects, where crops grew more vigorously due to better lighting, ventilation and population competition.

Canopy FPAR Measurements
The canola canopy FPAR measurements were performed after the UAV imag sition using the LP-80 PAR/LAI Ceptometer (AccuPAR, Decagon Devices Inc., P WA, USA). The LP-80 probe contained 80 independent sensors, which were evenly 1 cm apart. The photosensors measured the PAR in the 400 to 700 nm wa (µmol/(m 2 s)). During the FPAR measurement, the probe was placed along the ro tion within the canopy. Four positions and directions were considered to obtain hensive PAR data, i.e., the incident PAR, the top of canopy reflected PAR, the bo canopy transmission PAR and the bottom of soil reflection PAR. Each plot was r five times to calculate the average value and the canola canopy FPAR was then ca using Equation (1) [31]. A summary of the statistics of the canola FPAR acquired i given in Table 1. where Ic ↓ is the solar downward radiation at the canopy top, Ic ↑ is the solar radiati ted from the canopy, Iuc ↓ is the solar downward radiation under the canopy and I solar radiation reflected by the soil under the canopy.

Canopy FPAR Measurements
The canola canopy FPAR measurements were performed after the UAV image acquisition using the LP-80 PAR/LAI Ceptometer (AccuPAR, Decagon Devices Inc., Pullman, WA, USA). The LP-80 probe contained 80 independent sensors, which were evenly spaced 1 cm apart. The photosensors measured the PAR in the 400 to 700 nm waveband (µmol/(m 2 s)). During the FPAR measurement, the probe was placed along the row direction within the canopy. Four positions and directions were considered to obtain comprehensive PAR data, i.e., the incident PAR, the top of canopy reflected PAR, the bottom of canopy transmission PAR and the bottom of soil reflection PAR. Each plot was repeated five times to calculate the average value and the canola canopy FPAR was then calculated using Equation (1) [31]. A summary of the statistics of the canola FPAR acquired in situ is given in Table 1.
where I c↓ is the solar downward radiation at the canopy top, I c↑ is the solar radiation emitted from the canopy, I uc↓ is the solar downward radiation under the canopy and I uc↑ is the solar radiation reflected by the soil under the canopy.

PROSAIL RTM and Data Simulation
To investigate the optimal inversion method for estimating the FPAR based on the RTM, the crop canopy spectra were simulated using the PROSAIL model. The primary PROSAIL model was proposed by Baret et al. [32], who coupled the PROSPECT leaf optical properties model [33] with the SAIL canopy bidirectional reflectance model [34] to describe the variations in the canopy spectral reflectance as a function of the biochemistry and architecture parameters. At the leaf level, the PROSPECT model parameters included leaf chlorophyll a + b content (C ab ), carotenoid content (C ar ), brown pigment content (C brown ), equivalent water thickness (C w ), dry matter content (C m ) and leaf internal structure parameter (N). For the canopy level, the SAIL model parameters included the LAI, ALA, hot spot size parameters (hspot), SZA, observer zenith angle (OZA), relative azimuth angle (psi) and soil brightness parameters (psoil).
In general, the canopy structure parameters (e.g., LAI, SZA, ALA and psoil), as well as the leaf-level parameters (C ab and N), have been found to have an effect on the FPAR variation [12,35]. According to the prior information, the LAI, C ab , ALA and SZA were identified as particularly influential factors that affect the FPAR. To explore the effects of different parameters and ranges on the FPAR, two simulated datasets were generated. One was a general dataset parameter distribution according to the literature [36][37][38]. Based on prior information and the previous research results of Dong et al. [39], the specific dataset for canola was determined by reducing the value range and step size of the canopy structure parameters (e.g., C ab and LAI) and illumination parameter (SZA). The parameter ranges and specific meanings within the PROSAIL model for both datasets are shown in Table 2. Table 2. PROSAIL parameters range and specific meaning.

Model Parameters Typical Values
General Dataset 1 Range

Specific Dataset 2 Range
Step The FPAR could be derived by combining the law of conservation of energy with the PROSAIL model [40]. First, the leaf reflectance and transmittance were obtained using the PROSPECT part of the PROSAIL model. Second, the multiple scattering effect caused by the interaction between the canopy and background soil was partially accounted for by the SAIL component of the PROSAIL model. Finally, the canopy FPAR was solved by calculating both the direct and scattered PAR. The total FPAR, direct FPAR and diffuse FPAR were calculated using the following Equations (2)- (6): where α * s and α * d are the canopy absorptance for direct solar incident flux and hemispherical diffuse incident flux, respectively; α s and α d are the isolated canopy layer absorptance for direct solar incident flux and hemispherical diffuse incident flux respectively; τ ss and τ sd are the direct and bi-hemispherical transmittance for solar flux respectively; τ dd is the bi-hemispherical transmittance; r sd is the directional-hemispherical reflectance factor for solar incident flux and r dd is the bi-hemispherical reflectance factor; ρ sd and ρ dd are the directional-hemispherical reflectance for solar flux and bi-hemispherical reflectance above the canopy, respectively; ρ b dd is the bi-hemispherical reflectance at the bottom of the canopy; and E t dir and E t dif are the direct PAR and diffuse PAR, respectively.

Vegetation Indices (VIs)
Based on the simulated and measured canopy spectral data, a total of 29 VIs constructed using near-infrared and visible bands were selected ( Table 3). The selected VIs, including the NDVI, SR and difference vegetation index (DVI), as well as modified indices, were employed to investigate the correlation between the FPAR and VIs. In order to identify the VIs suitable for FPAR estimation using general and specific simulation datasets, the first step was to establish the relationship between the simulated VIs and the FPAR by employing curve-fitting models, such as linear, exponential, power and logarithm models. The second step was to determine the major VIs that exhibited a good correlation with the FPAR using estimated error metrics. Subsequently, the VIs were evaluated by calculating the rank sum and arranging them in ascending order accordingly. Table 3. Multispectral VIs for the FPAR estimation.

Sensitivity Analysis
To evaluate the influence of various leaf and canopy parameters on the variability of the FPAR, a global sensitivity analysis based on the extended Fourier amplitude sensitivity test (EFAST) algorithm was conducted. The EFAST is a variance-based global sensitivity analysis method that combines the computational efficiency of the classical FAST method with the global capability of the Sobol method. It is well-suited for nonlinear models, non-monotone models and interactions between parameters [68].
The EFAST analysis enabled the derivation of two indices: the first-and the total-order indices. The first-order index quantifies the individual contribution of each parameter to the overall variability, whereas the total-order index takes into account its interaction effect with other parameters. The 10,000 sample sets were generated based on the Monte Carlo method [69]. Subsequently, the FPAR was derived from the PROSAIL model and was then used to quantify the relative contribution of each model parameter.
In addition to the global sensitivity analysis, the local sensitivity of VIs to changes in the model parameters was further investigated. Appropriate VIs were not only sensitive to changes in the FPAR but also insensitive to other crop, soil and ambient factors. Hence, several sensitivity parameters determined beforehand using the EFAST method (e.g., LAI, SZA, ALA and psoil) were selected to evaluate the robustness of the VIs, while other insensitive parameters were kept fixed at the typical values (Table 2).

Inversion Modeling Algorithm
The regression relationship used in this study was established as an explicit parametric expression between the VIs and simulated FPAR from the PROSAIL model. The fitting functions could be linear or nonlinear (e.g., exponential or power function). Nonparametric algorithms, on the other hand, usually interpreted the strong nonlinearity of functional dependence between the biophysical/biochemical parameters and the spectral reflectance, such as ANN and SVR regression, which were more accurate and rational for the establishment of the inversion model. In this study, both the ANN and SVR algorithms were employed to establish a hybrid model with the selected VIs. For both the general and specific datasets, 80% of the samples were randomly extracted as the training (modeling) set, while the remaining samples were used for the model testing (validating) set. The process of establishing and validating hybrid inversion models based on the PROSAIL model and UAV images is presented in Figure 4. dependence between the biophysical/biochemical parameters and the spectral refle such as ANN and SVR regression, which were more accurate and rational for the lishment of the inversion model. In this study, both the ANN and SVR algorithm employed to establish a hybrid model with the selected VIs. For both the general an cific datasets, 80% of the samples were randomly extracted as the training (modelin while the remaining samples were used for the model testing (validating) set. The p of establishing and validating hybrid inversion models based on the PROSAIL mod UAV images is presented in Figure 4. The ANN algorithm mainly used the PROSAIL model as the training mod trained the ANN with the corresponding relationship between the canopy VIs a FPAR. After the training, the top of the canopy reflectance extracted from the UAV was used to evaluate the model performance. When applying the network model, put layer was the VIs and the output layer was the FPAR. Training an ANN neces determining the network type and structure (i.e., the number of hidden layers and per layer), and appropriately initializing the weights and regularization parame prevent overfitting. A two-layer feed-forward network with sigmoid hidden neuro linear output neurons using the Levenberg-Marquardt algorithm was optimize ANN weights were randomly initialized according to the Nguyen-Widrow metho a cross-validation procedure was implemented to mitigate overfitting [70].

The SVR Algorithm
The SVR is an empirical model based on a regression-type support vector m which can transform nonlinear regressions into linear ones by mapping the origin dimensional input space into a higher-dimensional feature space based on various functions [71]. To ensure its performance, the SVR required a training process that mined a set of model parameters by minimizing the generalization error. Prediction The ANN algorithm mainly used the PROSAIL model as the training model and trained the ANN with the corresponding relationship between the canopy VIs and the FPAR. After the training, the top of the canopy reflectance extracted from the UAV images was used to evaluate the model performance. When applying the network model, the input layer was the VIs and the output layer was the FPAR. Training an ANN necessitated determining the network type and structure (i.e., the number of hidden layers and nodes per layer), and appropriately initializing the weights and regularization parameters to prevent overfitting. A two-layer feed-forward network with sigmoid hidden neurons and linear output neurons using the Levenberg-Marquardt algorithm was optimized. The ANN weights were randomly initialized according to the Nguyen-Widrow method, and a cross-validation procedure was implemented to mitigate overfitting [70].

The SVR Algorithm
The SVR is an empirical model based on a regression-type support vector machine, which can transform nonlinear regressions into linear ones by mapping the original lowdimensional input space into a higher-dimensional feature space based on various kernel functions [71]. To ensure its performance, the SVR required a training process that determined a set of model parameters by minimizing the generalization error. Predictions were done in the SVR by using an optimal hyperplane of a Gaussian radial basis kernel to minimize the prediction error. The kernel parameters (e.g., penalty and loss function) were determined via a grid search based on the least mean square error.

Appropriate VIs for FPAR Estimation
To investigate the capability of VIs for interpreting the variations in the FPAR, the empirical regression models between different VIs and simulated FPAR based on two PROSAIL datasets are depicted in Table 4. The results indicate that the performances of different VIs for retrieving the FPAR varied across the two datasets, of which the parameter ranges were non-equivalent. Generally, the performance of the VIs could be summarized into four types of features. First, some VIs, like MCARVI, RDVI red-edge and MSAVI, presented good estimation performance for the general dataset but showed worse accuracies in the specific dataset, with a decline in R 2 and an increase in RMSE. Contrarily, certain VIs (e.g., NDVI, RVI, SR and mNDVI) displayed superior regression model performance only with a specific dataset. The third type of VIs demonstrated stable and consistent model capabilities for both general and specific datasets, including OSAVI, WDRVI, mSR and OSAVI red-edge . Lastly, the remaining VIs exhibited poor model properties compared with the others, indicating that these limited the effectiveness of capturing the FPAR under common conditions. Therefore, considering the modeling signatures and accuracy ranking for both datasets, the top three VIs of OSAVI, WDRVI and mSR, which presented robustness and adaptableness among the 29 VIs, were identified as effective variables for further retrieving the crop FPAR.

Global Sensitivity Analysis of the FPAR
A global sensitivity analysis was conducted using simulated canopy spectral data with the PROSAIL model and the EFAST algorithm to quantify the relative contributions of the leaf and canopy parameters to the variation in the FPAR ( Figure 5). The results consistently indicated that the three parameters of LAI, SZA and ALA produced the dominant influences on the variability of the FPAR among the nine parameters of interest, accounting for over 95% of total variance. C ab and psoil ranked fourth and fifth, respectively, in terms of their contributions. Although there was a slight difference in the total-order indices between the general and specific datasets, LAI remained the greatest contributor to the variations in the FPAR. For the general dataset, the ALA and SZA jointly dominated the variability of the FPAR, making up about 10% of the variance, except for the LAI. Meanwhile, as the parameters distribution narrowed (specific dataset), the sensi-tivity index of the LAI decreased from 0.79 to 0.67, and the sensitivity of ALA increased from 0.08 to 0.17. Although the influences of C ab and psoil on the FPAR variation were relatively weak, they were still much higher than those of the other four parameters (i.e., N, C m , hspot and C w ), of which the influences on the FPAR changes were less than one thousandth and could be considered negligible.
to the variations in the FPAR. For the general dataset, the ALA and SZA the variability of the FPAR, making up about 10% of the variance, e Meanwhile, as the parameters distribution narrowed (specific dataset dex of the LAI decreased from 0.79 to 0.67, and the sensitivity of ALA i to 0.17. Although the influences of Cab and psoil on the FPAR variat weak, they were still much higher than those of the other four para hspot and Cw), of which the influences on the FPAR changes were l sandth and could be considered negligible.

Local Sensitivity Analysis of the FPAR
In addition to the global sensitivity analysis, the variabilities of th functions of each parameter were investigated as well. Theoretically, not only have a significant relationship with the target variable (FPA resistance to changes in other variables. For the ALA parameter, the had a great influence on three VIs. The amount of canopy light interce ALA increased, resulting in a decline in the values of the VIs and F good linear regression between the VIs and FPAR could be observed less than 70°. Similar variabilities of OSAVI and mSR were achieved w decreasing ALA, but the WDRVI presented remarkable changes beyon (Figure 7a).

Local Sensitivity Analysis of the FPAR
In addition to the global sensitivity analysis, the variabilities of the VIs and FPAR as functions of each parameter were investigated as well. Theoretically, a robust VI should not only have a significant relationship with the target variable (FPAR) but also display resistance to changes in other variables. For the ALA parameter, the variations in ALA had a great influence on three VIs. The amount of canopy light interception decreased as ALA increased, resulting in a decline in the values of the VIs and FPAR (Figure 6a). A good linear regression between the VIs and FPAR could be observed when the ALA was less than 70 • . Similar variabilities of OSAVI and mSR were achieved with increasing and decreasing ALA, but the WDRVI presented remarkable changes beyond −100% and 100% (Figure 7a).
The leaf chlorophyll content was another important indicator for considering the sensitivity of the VIs. Figures 6b and 7b depict the variations in the VIs and FPAR with the C ab changes. In general, the VIs decreased with a drop in C ab due to an increase in the spectral reflectance of visible light. Compared with the mSR and WDRVI, the OSAVI presented excellent robustness for C ab changes. The relationships between OSAVI and FPAR were fully consistent in most cases, with variations only observed at a C ab of 10 µg cm −2 . The VIs of mSR and WDRVI showed good robustness with higher C ab . Meanwhile, the sensitivities of mSR and WDRVI to C ab variations were significantly increased when C ab was lower than 30 µg cm −2 . Substantial variations for mSR and WDRVI were noticed with a C ab of 10 µg cm −2 , implying that the estimation accuracy was exacerbated during the crop maturity period when the leaves turned yellow.
The effects of SZA changes on three VIs were illustrated in Figures 6c and 7c. The SZA had a minimal effect on the three VIs, except for slight changes at the lowest SZA (0 • ). Although the VIs presented a strong resistance to changes in SZA, the regressions between the VIs and FPAR were slightly different, with a high linear correlation with mSR and a nonlinear correlation with both OSAVI and WDRVI. These results indicated that the estimated FPAR suppressed the impact from SZA, which broke through the common time for capturing images at noon. The psoil was the important parameter in the PROSAIL, which could significantly influence the canopy spectral reflectance by controlling the soil brightness from dark (1) to bright (0). The variabilities in the VIs and FPAR in psoil are shown in Figures 6d and 7d. The VIs would be decreased with the increase in psoil. A smaller amplitude of regression for the variation in psoil was observed in OSAVI and the sensitivity of the psoil parameter was decreased under a higher canopy cover (FPAR > 0.8). Compared with OSAVI, WDRVI and mSR exhibited more variability and weaker resistance to changes in psoil, particularly when the FPAR was between 0.4 to 0.6.  The leaf chlorophyll content was another important indicator for considering the sensitivity of the VIs. Figures 6b and 7b depict the variations in the VIs and FPAR with the Cab changes. In general, the VIs decreased with a drop in Cab due to an increase in the spectral reflectance of visible light. Compared with the mSR and WDRVI, the OSAVI presented excellent robustness for Cab changes. The relationships between OSAVI and FPAR were fully consistent in most cases, with variations only observed at a Cab of 10 µg cm −2 . The VIs of mSR and WDRVI showed good robustness with higher Cab. Meanwhile, the sensitivities of mSR and WDRVI to Cab variations were significantly increased when Cab was lower than 30 µg cm −2 . Substantial variations for mSR and WDRVI were noticed with a Cab of 10 µg cm −2 , implying that the estimation accuracy was exacerbated during the crop maturity period when the leaves turned yellow.

Performance of the Inversion Model with the Validation Dataset
Although the robustness results of the VIs of OSAVI, WDRVI and mSR for interpreting the variability in FPAR were given in the previous section, the performance of the optimal estimation model was unstable when applied to different datasets. Therefore, in order to address this issue, hybrid inversion models based on the ANN and SVR were established to explore the optimal estimation of FPAR. The inversion performances of these models are displayed in Table 5. Consistent results were achieved showing that OSAVI had better predictive capability than the other VIs for all three inversion methods, with an outstanding R 2 and minimum RMSE. Among the three inversion methods, the ANN algorithm worked best, surpassing both the curve fitting and SVR methods. Overall, the ANN-OSAVI hybrid inversion model based on PROSAIL model data presented the optimal capability for estimating FPAR, with an R 2 of 0.822 and RMSE of 0.055 for general dataset 1 and an R 2 of 0.775 and RMSE of 0.059 for specific dataset 2.

Evaluation of the Optimal Inversion Model for Estimating Canola
Although the ANN-OSAVI-based hybrid model exhibited the optimal accuracy compared with WDRVI and mSR, the difference between the three hybrid models was minimal, from which it is challenging to determine the capability in estimating the canola FPAR due to potential discrepancies between the RTM model and real field conditions. Therefore, hybrid models of the ANN algorithm combined with three VIs were evaluated using the canola FPAR of each growing stage and the whole growing stages (Figure 8).
In general, the inversion accuracies of the hybrid model in the seedling stages were significantly higher than those in other stages. Similar trends were observed for the flowering and maturity stages. Specifically, for the seedling stage, ANN-OSVAI exhibited the highest accuracy, with R 2 and RMSE values of 0.78 and 0.043, respectively. ANN-WDRVI and ANN-mSR also performed well but had slightly lower R 2 values around 0.7 due to the discrete samples. Additionally, these models presented underestimated tendencies when the FPAR > 0.6 ( Figure 7a). With respect to the flowering stage, the ranges of FPAR were narrower and mainly concentrated between 0.5 and 0.7. The underestimation trend was greatly aggravated for all hybrid models, inducing increased estimated errors. ANN-OSAVI produced the closest fit to the 1:1 reference, followed by ANN-mSR and ANN-WDRVI. Overall, the ANN-OSAVI hybrid model presented the optimal capability for retrieving FPAR in canola growth periods, with the highest estimation accuracy of R 2 of 0.65 and an RMSE of 0.051. Based on the optimal inversion method, the estimation and inversion of the FPAR in different growth stages of winter canola were mapped, as shown in Figure 9. In general, the inversion accuracies of the hybrid model in the seedling stages were significantly higher than those in other stages. Similar trends were observed for the flowering and maturity stages. Specifically, for the seedling stage, ANN-OSVAI exhibited the  and an RMSE of 0.051. Based on the optimal inversion method, the estimation and inversion of the FPAR in different growth stages of winter canola were mapped, as shown in Figure 9.

Sensitivity Analysis of the FPAR and VIs
The purpose of the sensitivity analysis was to reveal the influences of the different leaf and canopy parameters on the variation in the crop FPAR and filter VIs with high resistance capability. For the EFAST, the range of PROSAIL model parameters also changed the contribution to the final FPAR. Previous studies showed that PROSAIL model parameters exhibited varying sensitivities across different spectral bands. For example, the red and near-infrared bands were sensitive to changes in Cab and LAI while being relatively insensitive to the soil matrix and atmospheric conditions [72]. ALA was sensitive to the near-infrared band, whereas the visible light band was more sensitive to SZA [73,74]. The results of two datasets with wide and narrow parameter ranges showed that the FPAR was mainly sensitive to the vegetation canopy structure parameters, including the LAI, ALA and SZA ( Figure 5). Both LAI and ALA determined the canopy absorption efficiency by adjusting the direct and diffuse directional transmittance and reflectance. SZA ultimately affected the FPAR simulation value by altering the total incoming PAR above the canopy [75]. Dong et al. [12] suggested that SZA, LAI, ALA, Cab and the canopy background reflectance had significant effects on the FPAR, while other parameters (hspot, Car, Cm) had lesser effects on the FPAR. The results of this study indicate that Cab and psoil also exerted certain effects on the FPAR.
The VIs are formed using a mathematical combination of different bands' reflectances to enhance the target signatures and eliminate other interferences [76,77]. Therefore, the relationships between the VIs and FPAR are comprehensively affected by various factors, such as atmospheric environment, vegetation type and remote sensing data quality [23,78,79]. For example, OSAVI eliminated effects from the soil background [80] and mSR might be less sensitive to canopy optical and geometric properties than NDVI [81]. WDRVI used NDVI near-infrared band reflectance for weighting calculations, reducing the weight of these bands under medium-and high-biomass conditions [82]. The sensitivity analysis of the VIs to the PROSAIL model parameters showed that the relationship between OSAVI and FPAR was less affected by psoil and Cab in comparison with mSR and WDRVI ( Figure 6). The findings regarding psoil agreed well with the proposal of Bannari and Staenz [83], but the results for Cab contradicted those reported by Liu et al. [84] and Bannari and Staenz [83]. This discrepancy can be attributed to the change in Cab content

Sensitivity Analysis of the FPAR and VIs
The purpose of the sensitivity analysis was to reveal the influences of the different leaf and canopy parameters on the variation in the crop FPAR and filter VIs with high resistance capability. For the EFAST, the range of PROSAIL model parameters also changed the contribution to the final FPAR. Previous studies showed that PROSAIL model parameters exhibited varying sensitivities across different spectral bands. For example, the red and nearinfrared bands were sensitive to changes in C ab and LAI while being relatively insensitive to the soil matrix and atmospheric conditions [72]. ALA was sensitive to the near-infrared band, whereas the visible light band was more sensitive to SZA [73,74]. The results of two datasets with wide and narrow parameter ranges showed that the FPAR was mainly sensitive to the vegetation canopy structure parameters, including the LAI, ALA and SZA ( Figure 5). Both LAI and ALA determined the canopy absorption efficiency by adjusting the direct and diffuse directional transmittance and reflectance. SZA ultimately affected the FPAR simulation value by altering the total incoming PAR above the canopy [75]. Dong et al. [12] suggested that SZA, LAI, ALA, C ab and the canopy background reflectance had significant effects on the FPAR, while other parameters (hspot, C ar , C m ) had lesser effects on the FPAR. The results of this study indicate that C ab and psoil also exerted certain effects on the FPAR.
The VIs are formed using a mathematical combination of different bands' reflectances to enhance the target signatures and eliminate other interferences [76,77]. Therefore, the relationships between the VIs and FPAR are comprehensively affected by various factors, such as atmospheric environment, vegetation type and remote sensing data quality [23,78,79]. For example, OSAVI eliminated effects from the soil background [80] and mSR might be less sensitive to canopy optical and geometric properties than NDVI [81]. WDRVI used NDVI near-infrared band reflectance for weighting calculations, reducing the weight of these bands under medium-and high-biomass conditions [82]. The sensitivity analysis of the VIs to the PROSAIL model parameters showed that the relationship between OSAVI and FPAR was less affected by psoil and C ab in comparison with mSR and WDRVI ( Figure 6). The findings regarding psoil agreed well with the proposal of Bannari and Staenz [83], but the results for C ab contradicted those reported by Liu et al. [84] and Bannari and Staenz [83]. This discrepancy can be attributed to the change in C ab content having a small effect on the FPAR at low LAI values. Moreover, this study did not observe the advantage of mSR being less sensitive to canopy optical and geometric properties (e.g., ALA and SZA). Xie et al. [85] also showed that mSR was susceptible to environmental impacts. WDRVI exhibited more fluctuations with C ab compared with OSAVI. In contrast, the anti-interference capability of WDRVI to parameter changes was the worst compared with OSAVI and mSR (Figure 7).

Comparison of Modeling Inversion Methods
Although the RTM is susceptible to the limitations of the inversion method and the complexity of the canopy radiation interaction process, it is considered a sounder mechanism for retrieving vegetation canopy signatures than the empirical model [86]. Among various inversion methods based on the RTM, the hybrid inversion method has the advantages of simplicity, accuracy and does not require initial values [40]. It has better potential for estimating crop biophysical and biochemical parameters when applied to remote sensing data [87]. In this study, the RTM combined with hybrid inversion was used to obtain the optimal method for estimating the FPAR. The results of the fitting regression showed that OSAVI, mSR and WDRVI exhibited higher correlations with the FPAR than other VIs and presented better performances in estimating the FPAR. Similar results were also reported by Leolini et al. [22], Dong et al. [88], and Viña and Gitelson [89]. Some classical VIs, such as SR, NDVI and RVI, performed well only with the specific dataset, and the overall rankings were slightly lower than with OSAVI, mSR and WDRVI. Therefore, although most VIs were shown to be effective in retrieving vegetation properties, their performances vary with different ranges [90].
It is necessary to determine the optimal VIs based on prior information. According to the results of Table 5, the ANN showed the best fitting effect, and the fitting accuracy of the SVR was higher than that of curve fitting but did not greatly improve the accuracy. Although curve fitting as an empirical model has high computational efficiency, it did not consider the complex transfer of solar radiation within the canopy, which may introduce large errors. The SVR was computationally efficient and optimized to be more accurate than other algorithms [91]. However, the SVR was more suitable for a small amount of sampled data and may become unpredictable when simulating a large amount of data, thereby weakening the predictive capability and ability to reflect the underlying patterns of the sample [92]. The ANN, on the other hand, has the greatest potential for solving nonlinear problems due to its accurate mapping capability [93]. The GEOV2 FPAR product uses ANN for FPAR estimation [94]. Therefore, when the ANN algorithm is combined with the VIs, it offers the advantages of easy interpretation, anti-interference and higher inversion accuracy.

Estimating the FPAR in Canola Growth Periods with Hybrid Models
This study demonstrated that the hybrid inversion models combining ANN and VIs were reliable for estimating the FPAR across the canola growth stages. However, comparing the FPAR estimated by the ANN-VIs hybrid models with the in situ measured FPAR, there was an underestimation trend at the flowering stage. Yellow flowers emitted red light radiation (yellow = green + red) for canopy horizontal signals, and this increased radiation caused an increase in the reflectance of the red band canopy [95,96]. The proportion of green and blue bands was positively correlated with the number of flowers per unit area [97]. This led to the estimated FPAR value being less than the measurement, inducing the lowest estimation for the canola flowering period. Furthermore, the stem gradually changed from green to withered yellow, while the pod gradually changed from green to yellow and gray, leading to errors in the estimation of the FPAR by the VIs at the maturity stage.
Traditional FPAR inversion studies mainly established correlations between the VIs and FPAR, such as SR and NDVI [98,99]. However, the relationships between the VIs and FPAR could be sensitive to change and unstable. Knyazikhin et al. [100] mentioned that the simple proportional assumption between the NDVI and FPAR was valid only when the canopy background was ideally black. There are also machine learning algorithms that can estimate the FPAR. Shi et al. [101] demonstrated the feasibility of estimating the maize FPAR using an ANN and stepwise multiple linear regression (SMLR) method (R 2 > 0.6). Although there were some errors in the FPAR estimation at the flowering and maturity stages, the hybrid model still exhibited acceptable performance during the seedling stage. This could provide a possibility for the accurate estimation of canola FPAR in the growth stages.

Conclusions and Recommendations
In this study, to identify a robust hybrid inversion model to improve the canola FPAR retrieving accuracy, the PROSAIL model was recompiled to simulate FPAR variations, and a total of 29 VIs and two machine learning algorithms were devoted to retrieving the FPAR. Three VIs of OSAVI, mSR and WDRVI had better capability in revealing the variations in the FPAR than the other VIs across both datasets. OSAVI, in particular, exhibited greater resistance to changes in the leaf and canopy parameters. The hybrid inversion model revealed that ANN-OSAVI exhibited the best performance for the canola FPAR retrieval, with whole R 2 and RMSE values of 0.65 and 0.051, respectively. In specific growth stages, the model achieved R 2 and RMSE values of 0.78 and 0.043 in the seedling stage, 0.63 and 0.047 in the flowering stage, and 0.61 and 0.050 in the mature stage, respectively.
This study highlights the potential of hybrid model methods in retrieving the canola FPAR during each growth stage. An accurately estimated FPAR could provide more veritable crop growth information for biomass prediction and vegetation productivity model (e.g., gross primary production (GPP)). For example, Clevers [102] estimated the time series of the FPAR of sugar beet and used it to calibrate the SUCROS model to achieve crop dry matter and yield predictions. Although these are encouraging results, further improvements are warranted. For instance, considering the internal structure of leaves, in addition to the woody part, there are other components that do not participate in photosynthesis but still absorb PAR, such as pigments, leaf veins and cell walls. The energy absorbed by these components accounts for 20% to 30% of the total leaf absorption [103]. Solar radiation absorbed by the non-photosynthetic part holds no practical significance for crop growth and development. Cheng et al. [2] investigated the impact of the LUE and FPAR changes on GPP estimation and proposed that the chlorophyll FPAR exhibited seasonal dynamics more similar to the flux tower than the canopy FPAR. Zhang et al. [104] studied the solar radiation absorption characteristics of trembling aspen forests to estimate the chlorophyll FPAR and canopy FPAR and suggested that the vegetation productivity had a strong correlation with the chlorophyll FPAR, but a weak correlation with the canopy FPAR. Therefore, future efforts should be devoted to developing and refining the FPAR from the canopy scale to the leaf scale to achieve higher accuracy in vegetation productivity estimations.
Author Contributions: Writing-original draft preparation, J.K.; methodology, J.K.; writing-review and editing, C.Z. and M.T.; data curation, Z.L. and R.L.; software, Z.X.; funding acquisition, C.Z. and S.F. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data that supported this study cannot be publicly shared due to ethical or privacy reasons and may be shared upon reasonable request to the corresponding author if appropriate.