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

Abstract

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 simplified 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 field 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 artificial 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² 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.

1. 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 field experiment were used to evaluate the inversion performance. The main objectives of this study were to (1) reveal the capability of different VIs in interpreting FPAR variation, (2) analyze the sensitivity of VIs to changes in leaf and canopy parameters, and (3) identify the robust hybrid estimation methods that consider inversion accuracy in different canola growth stages.

2. Materials and Methods

2.1. Study Area

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

2.2. 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⁻²), medium (MD 25 plants m⁻²) and high densities (HD 37.5 plants m⁻²) 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 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.

2.3. 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.

2.4. 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²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. Summary statistics of the canopy FPAR in three stages from 2020 to 2022.

Year	Growth Stages	N	Max	Min	Mean	SD
2021	Seedling	62	0.780	0.465	0.621	0.095
	Flowering	40	0.741	0.409	0.605	0.082
	Maturity	40	0.706	0.404	0.564	0.075
2022	Seedling	18	0.770	0.442	0.616	0.104
	Flowering	60	0.770	0.417	0.590	0.084
	Maturity	16	0.705	0.530	0.602	0.057

Note: The table reports the sample number (N), maximum value (max), minimum (min), mean (mean) and standard deviation (SD) of the observations.

.

$$\mathrm{FPAR}=\frac{({\mathrm{I}}_{\mathrm{c}\downarrow }-{\mathrm{I}}_{\mathrm{c}\uparrow }-({\mathrm{I}}_{\mathrm{uc}\downarrow }-{\mathrm{I}}_{\mathrm{uc}\uparrow }\left)\right)}{{\mathrm{I}}_{\mathrm{c}\downarrow }}$$

(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.

2.5. 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. PROSAIL parameters range and specific meaning.

Model	Parameters	Typical Values	General Dataset 1 Range	Specific Dataset 2 Range	Step
PROSPECT	Cab (μg cm−2)	40	10–80	30–60	0.5
	Car (μg cm−2)	8	5–20	5–20	10
	Cbrown (μg cm−2)	0.0	0–0.5	0–0.5	0.05
	Cw (cm)	0.01	0.005–0.05	0.005–0.05	0.005
	Cm (g cm−2)	0.009	0.001–0.015	0.001–0.015	0.001
	N	1.5	1–3	1–3	0.1
SAIL	LAI (m2 m−2)	1	0.25–7.5	1–5	0.5
	ALA (°)	30	10–80	30–60	10
	hspot (m m−1)	0.01	0–1	0–1	0.2
	SZA (°)	30	0–90	10–60	10
	OZA (°)	10	0–90	0–90	10
	psi (°)	0	0–90	0–90	10
	psoil	0	0–1	0–1	0.25

.

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):

$${\alpha }_{\mathrm{s}}^{*}={\alpha }_{\mathrm{s}}+\frac{{\mathrm{\tau }}_{\mathrm{ss}}{\mathrm{r}}_{\mathrm{sd}}+{\mathrm{\tau }}_{\mathrm{sd}}{\mathrm{r}}_{\mathrm{dd}}}{1-{\mathrm{r}}_{\mathrm{dd}}{\mathrm{\rho }}_{\mathrm{dd}}^{\mathrm{b}}}{\alpha }_{\mathrm{d}}$$

(2)

$${\alpha }_{\mathrm{d}}^{*}={\alpha }_{\mathrm{d}}+\frac{{\mathrm{\tau }}_{\mathrm{dd}}{\mathrm{r}}_{\mathrm{dd}}}{1-{\mathrm{r}}_{\mathrm{dd}}{\mathrm{\rho }}_{\mathrm{dd}}^{\mathrm{b}}}{\alpha }_{\mathrm{d}}$$

(3)

$${\alpha }_{\mathrm{s}}=1-{\mathrm{\rho }}_{\mathrm{sd}}-{\mathrm{\tau }}_{\mathrm{sd}}-{\mathrm{\tau }}_{\mathrm{ss}}$$

(4)

$${\alpha }_{\mathrm{d}}=1-{\mathrm{\rho }}_{\mathrm{dd}}-{\mathrm{\tau }}_{\mathrm{dd}}$$

(5)

$$\mathrm{FPAR}=\frac{\sum _{\mathsf{\lambda }=400}^{\mathsf{\lambda }=700}{(\alpha }_{\mathrm{s}}^{*}{\mathrm{E}}_{\mathrm{dir}}^{\mathrm{t}}+{\alpha }_{\mathrm{d}}^{*}{\mathrm{E}}_{\mathrm{dif}}^{\mathrm{t}})}{\sum _{\mathsf{\lambda }=400}^{\mathsf{\lambda }=700}({\mathrm{E}}_{\mathrm{dir}}^{\mathrm{t}}+{\mathrm{E}}_{\mathrm{dif}}^{\mathrm{t}})}={\mathrm{FPAR}}_{\mathrm{dir}}+{\mathrm{FPAR}}_{\mathrm{dif}}$$

(6)

where ${\alpha }_{\mathrm{s}}^{*}$ and ${\alpha }_{\mathrm{d}}^{*}$ are the canopy absorptance for direct solar incident flux and hemispherical diffuse incident flux, respectively; ${\alpha }_{\mathrm{s}}$ and ${\alpha }_{\mathrm{d}}$ are the isolated canopy layer absorptance for direct solar incident flux and hemispherical diffuse incident flux respectively; ${\mathrm{\tau }}_{\mathrm{ss}}$ and ${\mathrm{\tau }}_{\mathrm{sd}}$ are the direct and bi-hemispherical transmittance for solar flux respectively; ${\mathrm{\tau }}_{\mathrm{dd}}$ is the bi-hemispherical transmittance; ${\mathrm{r}}_{\mathrm{s}\mathrm{d}}$ is the directional-hemispherical reflectance factor for solar incident flux and ${\mathrm{r}}_{\mathrm{dd}}$ is the bi-hemispherical reflectance factor; ${\mathrm{\rho }}_{\mathrm{sd}}$ and ${\mathrm{\rho }}_{\mathrm{dd}}$ are the directional-hemispherical reflectance for solar flux and bi-hemispherical reflectance above the canopy, respectively; ${\mathrm{\rho }}_{\mathrm{dd}}^{\mathrm{b}}$ is the bi-hemispherical reflectance at the bottom of the canopy; and ${\mathrm{E}}_{\mathrm{dir}}^{\mathrm{t}}$ and ${\mathrm{E}}_{\mathrm{dif}}^{\mathrm{t}}$ are the direct PAR and diffuse PAR, respectively.

2.6. 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. Multispectral VIs for the FPAR estimation.

Vegetation Index	Formulation	References
Atmospherically resistant Vegetation index (ARVI)	$\frac{R_{840}-\left(2 \times R_{717}-R_{475}\right)}{R_{840}+\left(2 \times R_{717}-R_{475}\right)}$	[41]
Chlorophyll index green (CIgreen)	$\frac{R_{840}}{R_{560}}-1$	[42]
Chlorophyll index red-edge (CIred-edge)	$\frac{R_{840}}{R_{717}}-1$	[42]
Difference vegetation index (DVI)	$R_{840}-R_{664}$	[43]
Enhanced vegetation index (EVI)	$2.5 \times \frac{R_{840} - R_{668}}{R_{840} + {6 \times R}_{668} - 7.5 \times R_{475} + 1}$	[44]
Enhanced vegetation index 2 (EVI2)	$2.5 \times \frac{R_{840} - R_{668}}{1 + R_{840} + 2.4 \times R_{668}}$	[45]
Green normalized difference Vegetation index (GNDVI)	$\frac{R_{840} - R_{560}}{R_{840} + R_{560}}$	[46]
Green ratio vegetation index (RVIgreen)	$\frac{R_{840}}{R_{560}}$	[47]
Green-red vegetation index (GRVI)	$\frac{R_{560} - R_{668}}{R_{560} + R_{668}}$	[48]
Modified chlorophyll absorption ratio vegetation index (MCARVI)	$\left(\left(R_{840} - R_{717}\right) - 0.2 \times \left(R_{840} - R_{560}\right)\right) \times \left(\frac{R_{840}}{R_{717}}\right)$	[49]
Modified normalized difference vegetation index (mNDVI)	$\frac{R_{840} - R_{680}}{R_{840} + R_{680} - 2 \times R_{475}}$	[50]
modified normalized difference vegetation index red-edge (mNDVIred-edge)	$\frac{R_{840} - R_{717}}{R_{840} + R_{717} - 2 \times R_{475}}$	[51]
Modified simple ratio (mSR)	$\frac{\frac{R_{840}}{R_{668}} - 1}{\sqrt{\left(\frac{R_{840}}{R_{668}} + 1\right)}}$	[52]
Modified simple ratio red-edge (mSRred-edge)	$\frac{\frac{R_{840}}{R_{717}} - 1}{\sqrt{\left(\frac{R_{840}}{R_{717}} + 1\right)}}$	[53]
Modified triangular vegetation index 2 (MTVI2)	$\frac{\left(\right(1.5 \times (1.2 \times (R_{840} - R_{560}) - 2.5 \times (R_{668} - R_{560}\left)\right)}{\surd \left(\right({\left(2 \times R_{840} + 1\right)}^2 - (6 \times R_{840} - 5 \times \surd (R_{668}) - 0.5\left)\right)\left)\right)}$	[54]
Modified soil adjusted vegetation index (mSAVI)	$0.5 \times \left(2 \times R_{840} + 1 - {\left(2 \times R_{840} + 1\right)}^2 - 8 \times \sqrt{R_{840} - R_{668}}\right)$	[55]
Normalized difference vegetation index red-edge (NDVIred-edge)	$\frac{R_{840} - R_{717}}{R_{840} + R_{717}}$	[56]
Normalized difference vegetation index (NDVI)	$\frac{R_{840} - R_{668}}{R_{840} + R_{668}}$	[57]
Modified normalized difference water index (mNDWI)	$\frac{R_{560} - R_{840}}{R_{560} + R_{840}}$	[58]
Optimized soil-adjusted vegetation index (OSAVI)	$\frac{R_{840} - R_{668}}{R_{840} + R_{668} + 0.16}$	[59]
Optimized soil-adjusted vegetation index red-edge (OSAVIred-edge)	$\frac{R_{840} - R_{717}}{R_{840} + R_{717} + 0.16}$	[49]
Renormalized difference vegetation index (RDVI)	$\frac{R_{840} - R_{668}}{\sqrt{R_{840} + R_{668}}}$	[60]
Renormalized difference vegetation index red-edge (RDVIred-edge)	$\frac{R_{840} - R_{717}}{\sqrt{R_{840} + R_{717}}}$	[61]
Ratio vegetation index (RVI)	$\frac{R_{668}}{R_{840}}$	[62]
Soil-adjusted vegetation index (SAVI)	$\frac{1.5 \times R_{840} - R_{668}}{R_{840} + R_{668} + 0.5}$	[63]
Simple ratio (SR)	$\frac{R_{840}}{R_{668}}$	[64]
Simple ratio red-edge (SRred-edge)	$\frac{R_{717}}{R_{668}}$	[65]
Visible-band difference vegetation index (VDVI)	$\frac{{2 \times R}_{560}{ - R}_{668} - R_{475}}{2 \times R_{560} + R_{668} + R_{475}}$	[66]
Wide dynamic range vegetation index (WDRVI)	$\frac{{0.2 \times R}_{840} - R_{668}}{0.2 \times R_{840} + R_{668}}$	[67]

). 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.

2.7. 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).

2.8. 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.

2.8.1. The ANN Algorithm

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].

2.8.2. 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 low-dimensional 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.

3. Results

3.1. 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. Regression equation derived from simulated general and specific datasets.

VI	General Dataset 1				Specific Dataset 2				Total Rank
	Regression Equation	R2	RMSE	Rank	Regression Equation	R2	RMSE	Rank
OSAVI	$y=1.217x^{0.807}$	0.817	0.054	1	$y=1.164x^{0.634}$	0.771	0.100	5	6
WDRVI	$y=0.405x+0.650$	0.766	0.061	4	$y=0.463x+0.648$	0.823	0.088	2	6
mSR	$y=0.556x^{0.335}$	0.766	0.061	6	$y=0.559x^{0.344}$	0.818	0.089	4	10
MTVI2	$y=0.983x^{0.412}$	0.776	0.060	2	$y=0.977x^{0.335}$	0.712	0.112	9	11
RVI	$y=0.4548x^{-0.227}$	0.753	0.063	10	$y={0.253x}^{-0.391}$	0.821	0.089	3	13
NDVI	$y=0.910x+0.055$	0.744	0.064	12	$y=0.893x+0.108$	0.840	0.084	1	13
OSAVIred-edge	$y=1.442x^{0.634}$	0.764	0.057	8	$y=1.315x^{0.475}$	0.726	0.109	8	16
MCARVI	$y=0.883x^{0.237}$	0.769	0.061	3	$y=0.895x^{0.195}$	0.689	0.117	13	16
SR	$y=0.455x^{0.227}$	0.753	0.063	11	$y=0.431x^{0.247}$	0.756	0.103	6	17
RDVIred-edge	$y=1.332x^{0.565}$	0.766	0.061	5	$y=1.259x^{0.458}$	0.696	0.115	12	17
RDVI	$y=1.187x^{0.646}$	0.754	0.063	9	$y=1.133x^{0.517}$	0.707	0.113	10	19
MSAVI	$y=1.029x^{0.488}$	0.755	0.063	7	$y=1.021x^{0.403}$	0.665	0.121	15	22
mNDVI	$y=0.908x^{1.580}$	0.683	0.071	17	$y=0.934x^{1.208}$	0.747	0.105	7	24
NDVIred-edge	$y=1.189x^{0.637}$	0.691	0.070	16	$y=1.153x^{0.505}$	0.703	0.114	11	27
SAVI	$y=1.159x^{0.588}$	0.729	0.066	13	$y=1.097x^{0.439}$	0.656	0.123	19	32
EVI2	$y=1.038x^{0.511}$	0.728	0.066	14	$y=1.018x^{0.412}$	0.659	0.122	18	32
GNDVI	$y=1.175x-0.037$	0.671	0.073	18	$y=2.289x-0.898$	0.673	0.120	14	32
EVI	$y=1.014x^{0.508}$	0.710	0.068	15	$y=0.997x^{0.402}$	0.630	0.127	20	35
mNDWI	$y=-1.175x-0.037$	0.671	0.073	19	$y=-1.077x+0.067$	0.660	0.122	17	36
mSRred-edge	$y=0.758x^{0.403}$	0.665	0.073	20	$y=0.793x^{0.340}$	0.664	0.121	16	36
CIred-edge705	$y=0.596x^{0.285}$	0.633	0.077	21	$y=0.633x^{0.245}$	0.594	0.133	22	43
CIgreen	$y=0.496x^{0.287}$	0.604	0.080	23	$y=0.529x^{0.256}$	0.563	0.138	23	46
SRred-edge	$y=0.476x^{0.415}$	0.610	0.079	22	$y=0.506x^{0.373}$	0.548	0.140	25	47
DVI	$y=1.187x^{0.394}$	0.599	0.080	24	$y=1.133x^{0.327}$	0.549	0.140	24	48
GRVI	$y=0.522x+0.635$	0.517	0.088	28	$y=0.666x+0.595$	0.600	0.132	21	49
mNDVIred-edge	$y=1.147x^{0.731}$	0.538	0.086	26	$y=1.105x^{0.561}$	0.543	0.141	26	52
GRVI	$y=0.427x^{0.340}$	0.581	0.082	25	$y=0.457x^{0.304}$	0.520	0.145	28	53
ARVI	$y=0.824x+0.597$	0.518	0.088	27	$y=0.708x+0.635$	0.415	0.160	29	56
VDVI	$y=0.581x+0.577$	0.472	0.092	29	$y=0.717x+0.529$	0.543	0.141	27	56

. 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² 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.

3.2. Sensitivity Analysis

3.2.1. 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 sensitivity 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.

3.2.2. 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. Figure 6b and Figure 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⁻². 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⁻². Substantial variations for mSR and WDRVI were noticed with a C_ab of 10 μg cm⁻², 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 Figure 6c and Figure 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 Figure 6d and Figure 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.

3.3. 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. Comparison of the hybrid inversion models for estimating FPAR.

Dataset	Inversion Method	VIs	Model Performance
			R2	RMSE
General dataset 1	Curve fitting	OSAVI	0.801	0.077
		mSR	0.703	0.122
		WDRVI	0.750	0.080
	ANN	OSAVI	0.822	0.055
		mSR	0.772	0.069
		WDRVI	0.769	0.062
	SVR	OSAVI	0.817	0.067
		mSR	0.741	0.088
		WDRVI	0.768	0.072
Specific dataset 2	Curve fitting	OSAVI	0.740	0.073
		mSR	0.680	0.120
		WDRVI	0.735	0.080
	ANN	OSAVI	0.775	0.059
		mSR	0.750	0.079
		WDRVI	0.752	0.066
	SVR	OSAVI	0.749	0.066
		mSR	0.722	0.083
		WDRVI	0.737	0.069

. Consistent results were achieved showing that OSAVI had better predictive capability than the other VIs for all three inversion methods, with an outstanding R² 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² of 0.822 and RMSE of 0.055 for general dataset 1 and an R² of 0.775 and RMSE of 0.059 for specific dataset 2.

3.4. 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² and RMSE values of 0.78 and 0.043, respectively. ANN-WDRVI and ANN-mSR also performed well but had slightly lower R² 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² 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.

4. Discussion

4.1. 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 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).

4.2. 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.

4.3. 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² > 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.

5. 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² and RMSE values of 0.65 and 0.051, respectively. In specific growth stages, the model achieved R² 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.