Evaluating Empirical Regression, Machine Learning, and Radiative Transfer Modelling for Estimating Vegetation Chlorophyll Content Using Bi-Seasonal Hyperspectral Images

Abstract

Different types of methods have been developed to retrieve vegetation attributes from remote sensing data, including conventional empirical regressions (i.e., linear regression (LR)), advanced empirical regressions (e.g., multivariable linear regression (MLR), partial least square regression (PLSR)), machine learning (e.g., random forest regression (RFR), decision tree regression (DTR)), and radiative transfer modelling (RTM, e.g., PROSAIL). Given that each algorithm has its own strengths and weaknesses, it is essential to compare them and evaluate their effectiveness. Previous studies have mainly used single-date multispectral imagery or ground-based hyperspectral reflectance data for evaluating the models, while multi-seasonal hyperspectral images have been rarely used. Extensive spectral and spatial information in hyperspectral images, as well as temporal variations of landscapes, potentially influence the model performance. In this research, LR, PLSR, RFR, and PROSAIL, representing different types of methods, were evaluated for estimating vegetation chlorophyll content from bi-seasonal hyperspectral images (i.e., a middle- and a late-growing season image, respectively). Results show that the PLSR and RFR generally performed better than LR and PROSAIL. RFR achieved the highest accuracy for both images. This research provides insights on the effectiveness of different models for estimating vegetation chlorophyll content using hyperspectral images, aiming to support future vegetation monitoring research.

1. Introduction

Vegetation biochemical and biophysical properties, such as chlorophyll content and leaf area index, are essential vegetation characteristics that influence plant physiological status, vegetation productivity, and ecosystem health [1,2]. As an example, chlorophyll controls the amount of solar radiation and CO₂ that a plant uses for photosynthesis, thus highly influencing vegetation photosynthetic capacity and productivity [3]. Investigation of vegetation properties is critical for understanding vegetation growth condition and supporting resource management [4]. Field measurements of vegetation properties are costly, labor-intensive, and limited to small areas [5]. In contrast, remote sensing is a low-cost and efficient tool for estimating vegetation properties and investigating their spatio-temporal variations over large areas [6,7]. Different approaches have been applied for estimating vegetation properties from remote sensing data, including conventional empirical regressions (i.e., linear regression (LR)), advanced empirical regressions (e.g., multivariable linear regression (MLR), partial least square regression (PLSR), principal component regression (PCR)), machine learning (e.g., random forest regression (RFR), decision tree regression (DTR), artificial neural networks (ANN)), and radiative transfer modelling (RTM, e.g., PROSPECT, PROSAIL) [5,8,9,10,11,12,13]. These approaches were developed using different theories and thus are expected to perform differently when estimating vegetation properties.

Previous research has attempted to compare various models for estimating vegetation properties. For instance, Powell et al. [14] compared RFR, reduced major axis regression, and gradient nearest neighbour imputation for estimating forest biomass using Landsat imagery. Darvishzadeh et al. [15] investigated PLSR and PROSAIL for retrieving grassland leaf area index (LAI) from HyMap data. Siegmann and Jarmer [16] evaluated performance of PLSR, support vector regression (SVR), and RFR for estimating wheat LAI using field-measured hyperspectral reflectance. Wang et al. [17] compared RFR, SVR, and ANN for investigating wheat biomass using multispectral HJ-CCD imagery. Reddy et al. [18] applied PLSR and RFR for modelling forest structural attributes (e.g., volume, tree height) using airborne multispectral imagery. Xing et al. [19] tested performance of PLSR, SVM, RFR, and ridge regression for estimating nutrition components (e.g., protein and sugar) in forage crops using a few field-based multispectral sensors and active sensors (e.g., laser). Yue et al. [13] examined ANN, MLR, DTR, PLSR, RFR, SVM, boosted binary regression tree (BBRT), and PCR for estimating wheat biomass using field-measured hyperspectral reflectance. In summary, most of these studies have compared advanced empirical regressions (e.g., PLSR, MLR) with machine learning methods (e.g., RFR, ANN), while less studies have compared these with two other commonly used approaches: conventional empirical regressions and radiative transfer modelling. The conventional empirical regression (e.g., LR using vegetation index) is the most widely used approach for retrieving vegetation information from remote sensing data owing to its easy application and fast computing. The advanced empirical regressions (e.g., PLSR) and machine learning approaches (e.g., RFR) are expected to outperform the conventional empirical regressions (e.g., LR) since they can use information from multiple predictor variables [20], however, it is unknown how much higher accuracy they can achieve. In many real-world practices, e.g., estimating crop properties in precision farming, the LR with acceptable accuracy is still preferred owing to its operational and computing advantages. Therefore, it is valuable to explore strategies for building well-performing LRs and comparing their accuracy with that of more advanced models. RTM is another approach used for investigating vegetation properties using remote sensing data. It describes the interaction of electromagnetic radiation with plants (reflection, transmission, absorption) using physical laws and simulates vegetation spectra at different wavelength ranges using vegetation biophysical and biochemical properties [21,22]. The RTM is a less often used approach due to model complexity and computing challenges. However, the RTM is physically-based, not sensor-, site-, or season-specific, and thus it is a more transferable approach than the empirical-based methods [5]. Therefore, it is valuable to compare RTM to other approaches and evaluate their accuracy and effectiveness.

Previous studies comparing different approaches have mainly utilized multispectral images or ground based hyperspectral reflectance data, while few studies have applied hyperspectral images to evaluate the performance of different models. The hyperspectral imagery records the spectral signal of ground features in hundreds of narrow spectral bands over a large area. The huge amount of spectral and spatial (i.e., image textural) information may influence the performance of different models, which has not been fully explored in previous studies. Moreover, most of previous studies have utilized single-date images for the comparison of different approaches, while less have used multi-temporal data. The phenological variations of vegetation (e.g., changes from homogeneous green canopies in the middle growing season to heterogeneous canopies with mixed green and senescent vegetation in the late growing season) also potentially influence the effectiveness of different models. It is thus critical to investigate if the models perform differently for images acquired from different seasons. In addition, most of previous model comparison studies attempted to investigate vegetation biophysical properties [14,15,16,17,18], such as vegetation structure and biomass, while few have focused on estimating vegetation biochemical properties (e.g., vegetation chlorophyll content). Different spectral features need to be utilized for investigating vegetation biochemical properties and model performance may vary depending on which spectral features are used. For instance, spectral signals in blue, red, and NIR ranges are critical for estimating canopy chlorophyll content (e.g., signals in blue and red capturing chlorophyll absorption features and signals in NIR reflecting vegetation density information). Models capable of selecting and utilizing spectral variables in these ranges can potentially achieve higher estimation accuracy. In this study, bi-seasonal high-spatial resolution hyperspectral images were acquired by a helicopter-based imaging system and utilized for estimating vegetation canopy chlorophyll content with different models. Specifically, LRs, PLSR, RFR, and a modified PROSAIL, which represent conventional empirical regression, advanced empirical regression, machine learning, and RTM, respectively, were selected and compared. Different types of information, such as vegetation indices (VIs), reflectance of each band, principal components (PCs), and textural variables were extracted from the hyperspectral images and applied into different models as predictor variables.

LR is the most widely used approach for estimating vegetation properties from remote sensing data [4]. One common practice is to establish a LR between targeted vegetation properties (e.g., chlorophyll content) and remote sensing variables (e.g., VIs) using sampling data, and then applying this LR for predictive purposes. A large number of broadband or narrowband VIs have been developed in previous studies to retrieve vegetation properties, such as the widely used Normalized Difference Vegetation Index (NDVI), Soil-Adjusted Vegetation Index (SAVI), and Atmospherically Resistant Vegetation Index (ARVI) [23,24,25,26]. Previous studies have also compared these VIs and evaluated their performance for estimating various vegetation properties (e.g., chlorophyll content, biomass) in different ecosystems (e.g., forests, grasslands, and croplands) [4,5,17,27,28,29,30]. In this research, 29 VIs that performed well in previous studies for estimating vegetation canopy chlorophyll content were selected and further evaluated [28,30]. In addition to VIs, previous studies indicated that image PCs can provide important information for investigating vegetation properties [31,32]. Therefore, PCs were also tested as predictor variables in LR for estimating vegetation chlorophyll content.

To estimate vegetation properties, LR is relatively easy to build in comparison with more advanced regression models and LR generally performs well with appropriate VIs and training samples (e.g., selecting VIs from previous studies, collecting samples in a wide data range). However, there are several problems associated with LR for the retrieval of vegetation properties, such as saturation of some VIs (e.g., NDVI for dense green canopies), being highly affected by sample data (e.g., extreme values), and VIs influenced by environment factors (e.g., soil background, view geometry, atmospheric condition) [33,34,35]. Advanced empirical regressions, such as PLSR, have been developed and utilized to deal with these problems and improve the stability of regression models [13,16]. PLSR is a bi-linear calibration method using data compression technique that reduces the number of collinear predictor variables to a few non-correlated latent variables [18,35]. These latent variables are generated to maximize information content from the original variables and to optimize their explanatory power for predictions [36]. PLSR is considered as an efficient method for feature extraction and dimension reduction without losing much information [18,37]. It combines features of principal component regression and stepwise multivariable regression, and is capable of tackling some common problems in regression models, such as collinearity, data noise, and over-fitting [13,37,38,39]. PLSR has been used in previous studies for estimating vegetation LAI, biomass, or nitrogen content [8,35,39,40]. In this study, PLSR was tested to estimate canopy chlorophyll content.

Random Forest is an ensemble modelling technique that builds a forest with a large number of regression trees using bootstrapped training data [41,42]. It is a robust and widely used prediction model for regression and classification [12,43]. It has many advantages over other regression models, such as being able to handle thousands of input variables and running efficiently on a large dataset, not being sensitive to noise or over-fitting, handling the collinearity problem, requiring a minimal number of parameters, and making no distribution assumptions about the predictor or response variables [13,17,44]. One essential feature of random forest is that the importance of predictor variables can be evaluated in the model, which is essential for researchers to understand the contributions of different variables to the regression or classification model [44]. This variable importance is assessed based on how much the prediction accuracy decreases when a variable is permuted while the others are left unchanged [45,46]. Random forest has been widely applied in the remote sensing field for classification [43,47], however, it has been less used for regression to estimate vegetation properties [17,33,34]. A few existing studies include estimating biomass [17], vegetation structure [44], and water content [48]. In this study, RFR is examined to investigate its performance at estimating vegetation canopy chlorophyll content.

PROSAIL is one of the most-widely applied canopy-level radiative transfer models, describing the spectral and directional variation of canopy reflectance in the solar domain [11]. It integrates the leaf-level PROSPECT model and canopy-level SAIL model. Specifically, it passes the output leaf reflectance and transmittance of the PROSPECT model into the SAIL model to simulate the whole spectro-directional canopy reflectance. Therefore, it links the canopy reflectance with leaf biochemical contents (e.g., leaf pigments, water, and dry matter), canopy architecture (e.g., LAI, leaf angle distribution, and relative leaf size), soil properties, illumination conditions, and viewing geometry [49]. The PROSAIL model has been applied successfully in different homogeneous ecosystems (e.g., forests and crops) for estimating vegetation biophysical and biochemical properties (e.g., LAI, leaf chlorophyll content) [11,50,51,52]. However, it was found to struggle with simulating spectra of heterogeneous canopies with mixed green and senescent vegetation [22,53], likely because the leaf-level PROSPECT model cannot sufficiently simulate spectra of senescent leaves [54]. A modified PROSPECT that is capable of tackling this problem was developed and applied in this study to generate a modified PROSAIL [54]. This modified PROSAIL is compared to regression models for estimating vegetation canopy chlorophyll content.

Canopy heterogeneity (i.e., mixed green and senescent vegetation) is an essential vegetation biophysical characteristic of an ecosystem, especially for heterogeneous grasslands or wetlands. Such heterogeneity may influence the retrieval of vegetation properties using remote sensing. As mentioned previously, high heterogeneity was found to bring considerable challenges to the PROSAIL model for simulating vegetation spectra in a grassland area [22,53]. In addition, it also affected the accuracy of empirical models (e.g., PLS) for estimating vegetation properties [55]. The influence of vegetation heterogeneity on the performance of different models (RTM-based or empirical models) for the estimation of vegetation properties have rarely been explored in previous research. Many studies estimating vegetation properties using different models have focused on homogeneous ecosystems, such as crops and forests, while less focused on heterogeneous ecosystems, such as grasslands and wetlands [14,56,57,58]. In this study, bi-seasonal hyperspectral images were acquired in a grassland area, which capture vegetation information with different level of heterogeneities. Specifically, one image was acquired in the middle growing season (i.e., June) that featured by homogeneous green canopies and the other image was acquired in the late growing season (i.e., August) that exhibited high heterogeneity (e.g., a mixture of green and senescent vegetation). These two images were used for evaluating effects of heterogeneity on the performance of models.

This study compared the performance of LR, PLSR, RFR, and a modified PROSAIL model for retrieving vegetation chlorophyll content from bi-seasonal hyperspectral images in a heterogeneous grassland area. Different image features, including VIs, spectral reflectance, PCs, and image textural variables were extracted from the images, aiming to acquire different types of information from the images and optimize the performance of different models. Contributions of different types of variables to different models (e.g., importance values of these variables) were also evaluated. Different optimization techniques were used in PLSR and RFR models for improving their performance. A modified PROSAIL model, which can accurately simulate reflectance of heterogeneous canopies, was adopted in this study and compared to other models. Factors influencing model performance were also discussed in order to provide insights with other researchers on the selection and optimization of different models.

2. Material and Methods

2.1. Study Area

This research was performed at Koffler Scientific Reserve (KSR), located in Southern Ontario, Canada (Figure 1). Grasslands in KSR mainly consist of temperate tall grasses with a height ranging from 10 to 50 cm during the growing season (May to September) [59]. Dominant grass species in this area are Awnless brome (Bromus inermis), Fescue (Festuca rubra L.), Goldenrod (Solidago canadensis L.), and Milkweed (Asclepias L.). The climate type in this area is temperate continental, with a mean temperature ranging from −10 °C (February) to 30 °C (July), and precipitation ranging from 20 mm (March) to 100 mm (July) [60]. The soil in this area is typically grey-brown podzolic soil.

2.2. Hyperspectral Imagery Collection and Field Survey

Hyperspectral images were collected using a Micro-HyperSpec (Headwall Photonics Inc., Boston, MA, USA) that was mounted on a manned helicopter (Figure 2). This sensor is capable of collecting spectral signals in more than 300 bands ranging from 400 to 1000 nm. Two flight missions were conducted around noon on 14 June (middle growing season) and 26 August, 2016 (late growing season), aiming to capture the homogeneous and heterogeneous features of this grassland ecosystem, respectively. The weather conditions on both days were sunny and clear. The flights were operated at a height of 250 m and acquired imagery with a spatial resolution of 30 cm. The acquired images were radiometrically and geometrically corrected using SpectralView that was provided by the sensor manufacturer. The images were further atmospherically corrected in ENVI (Exelis Visual Information Solutions, Boulder, CO, USA) using an empirical line method [61,62]. Lastly, the images were resampled to 301 bands with a 2 nm interval from 400 to 1000 nm.

Field surveys were performed simultaneously with the flight missions. A total of 29 study sites were pre-selected in the study area for the field survey in June. The selection of sites was dependent on vegetation growing conditions, species composition, and topographic conditions. An extra five study sites were added for the field survey in August due to the increased variety of vegetation canopy features. The location of each site was obtained using a highly accurate Trimble GeoExplorer GPS (Trimble Navigation Limited, Sunnyvale, CA, USA). Field data collected at the study sites included spectral reflectance, LAI, vegetation height, species composition, and canopy photos. The spectral reflectance data were collected using an ASD spectroradiometer FieldSpec 3 in the range of 350–2500 nm (Analytical Spectral Devices Inc., Boulder, CO, USA). LAI was measured using an AccuPAR LP-80 ceptometer (Decagon Devices, Inc., Pullman, Washington, WA, USA). Leaf samples for each site were also collected and transported to a laboratory for measuring leaf reflectance and chlorophyll content. Leaf reflectance was acquired using Plant Probe (Analytical Spectral Devices, Inc., Boulder, CO, USA) that was connected to the ASD FieldSpec spectroradiometer, following the protocol provided by [63] and [64]. The measured reflectance data were used for evaluating the simulation results from radiative transfer models. Leaf chlorophyll was extracted using N, N-dimethylformamide (DMF) following a protocol proposed by [65]. Canopy chlorophyll content was then calculated using leaf chlorophyll content times LAI of green vegetation (i.e., green LAI) [2,5,66,67,68]. For the survey in June, the green LAI equals to the field-measured canopy LAI since the canopies were homogeneous green vegetation. For the survey in August, the canopies were heterogeneous with mixed green and senescent vegetation. The green LAI was calculated using field-measured canopy LAI times the percentage of green vegetation in this canopy. This percentage (e.g., 20%, 55%, or 100%) was manually estimated in the field using a quadrat (0.5 × 0.5 m). Statistics of measured leaf chlorophyll content and green LAI is shown in Figure 3. Correlations between canopy chlorophyll content and green LAI and correlations between canopy chlorophyll and leaf chlorophyll were calculated to evaluate the influence of green LAI and leaf chlorophyll on the variations of canopy chlorophyll.

2.3. Methods and Model Parameter Settings

LR, PLSR, RFR, and a modified PROSAIL were selected for comparison to evaluate their performance of estimating canopy chlorophyll content. Details of each model are described in below sections. A leave-one-out cross validation (LOOCV) was used for validating the models [13,33] and the coefficient of determination (R²), together with root mean square error (RMSE), were calculated and applied to describe performance of different models [13].

2.3.1. Linear Regression

LRs were built using different predictor variables to estimate canopy chlorophyll content. In this study, 29 narrowband VIs (

Table 1. Narrowband vegetation indices selected in this research. (R in the formula indicates reflectance, values indicate wavelength in nm).

Index	Full Name	Formula	References
BGI	Blue/Green Pigment Index	$R_{450}/R_{550}$	[70]
DVI	Difference Vegetation Index	$R_{800}-R_{680}$	[71]
GI	Greenness Index	$R_{554}/R_{677}$	[70]
GVI	Greenness Vegetation Index	$(R_{682}-R_{553})/\left(R_{682}+R_{553}\right)$	[72]
MCARI1	Modified Chlorophyll Absorption Ratio Index 1	$\left[\left(R_{750}-R_{705}\right)-0.2\left(R_{750}-R_{550}\right)\right]\left(R_{750}/R_{705}\right)$	[73]
MCARI2	Modified Chlorophyll Absorption Ratio Index 2	$\frac{1.5\left[1.2\left(R_{800}-R_{550}\right)-2.5\left(R_{670}-R_{550}\right)\right]}{\sqrt{{\left(2R_{800}+1\right)}^2-\left(6R_{800}-5\sqrt{R_{670}}\right)}-0.5}$	[74]
MRENDVI	Modified Red Edge Normalized Difference Vegetation Index	$\left(R_{750}-R_{705}\right)/\left(R_{750}+R_{705}-2R_{445}\right)$	[75]
MSAVI	Modified Soil Adjusted Vegetation Index	$0.5\left[2R_{800}+1-\sqrt{{\left(2R_{800}+1\right)}^2-8\left(R_{800}-R_{670}\right)}\right]$	[25]
mSR1	Modified Simple Ratio 1	$\frac{\frac{R_{800}}{R_{670}}-1}{\sqrt{\frac{R_{800}}{R_{670}}}+1}$	[76]
mSR2	Modified Simple Ratio 2	$\frac{\left(\frac{R_{750}}{R_{705}}\right)-1}{\sqrt{\left(\frac{R_{750}}{R_{705}}\right)+1}}$	[73]
mSR3	Modified Simple Ratio 3	$\left(R_{750}-R_{445}\right)/\left(R_{705}-R_{445}\right)$	[75]
MTCI	MERIS Terrestrial Chlorophyll Index	$\left(R_{754}-R_{709}\right)/\left(R_{709}-R_{681}\right)$	[77]
MTVI	Modified Triangular Vegetation Index	$\frac{1.5\left[1.2\left(R_{800}-R_{550}\right)-2.5\left(R_{670}-R_{550}\right)\right]}{\sqrt{{\left(2R_{800}+1\right)}^2-6R_{800}+5{\left(R_{670}\right)}^{0.5}-0.5}}$	[74]
NDRE	Normalized Difference Red-edge index	$\left(R_{790}-R_{720}\right)/\left(R_{790}+R_{720}\right)$	[78]
NDVI	Normalized Difference vegetation index	$\left(R_{682}-R_{553}\right)/\left(R_{682}+R_{553}\right)$	[72]
OSAVI1	Optimized Soil Adjusted Vegetation Index 1	$1.16\left(R_{800}-R_{670}\right)/\left(R_{800}+R_{670}+0.16\right)$	[79]
OSAVI2	Optimized Soil Adjusted Vegetation Index 2	$1.16\left(R_{750}-R_{705}\right)/\left(R_{750}+R_{705}+0.16\right)$	[73]
PPR	Plant Pigment Ratio	$\left(R_{550}-R_{450}\right)/\left(R_{550}+R_{450}\right)$	[80]
PRI	Photochemical Reflectance Index	$\left(R_{570}-R_{539}\right)/\left(R_{570}+R_{539}\right)$	[81]
RDVI	Renormalized Difference vegetation index	$\frac{R_{800}-R_{670}}{\sqrt{R_{800}+R_{670}}}$	[82]
RENDVI	Red Edge Normalized Difference Vegetation Index	$\left(R_{750}-R_{705}\right)/\left(R_{750}+R_{705}\right)$	[75,83]
REPI	Red Edge Position Index	$700+40\left({\mathsf{\rho }}_{\mathrm{REP}}-R_{700}\right)/\left(R_{740}-R_{700}\right)$ ${\mathsf{\rho }}_{\mathrm{REP}}=0.5\left(R_{670}+R_{780}\right)$	[84]
RRVI	Reciprocal Reflectance-based Vegetation Index	$R_{750-800}/R_{695-740}-1$	[85]
RGI	Red/Green Index	$R_{690}/R_{550}$	[70]
SPVI	Spectral Polygon Vegetation Index	$0.4\left[3.7\left(R_{800}-R_{670}\right)-1.2\left|R_{530}-R_{670}\right|\right]$	[86]
SR	Simple Ratio	$R_{800}/R_{675}$	[71]
TSAVI	Transformed Soil Adjusted Vegetation Index	$\frac{\alpha \left(R_{875}-\alpha R_{680}-\beta \right)}{\left[R_{680}+\alpha \left(R_{875-\beta }\right)+0.08\left(1+{\alpha }^2\right)\right]}$ $\alpha =1.062\beta =0.022$	[26]
VREI	Vogelmann Red Edge Index	$\left(R_{734}-R_{747}\right)/\left(R_{715}+R_{726}\right)$	[87]
ZM	Zarco and Miller	$R_{750}/R_{710}$	[51]

) that performed well in previous studies for estimating canopy chlorophyll content were used as predictor variables [28,30]. In addition, a preliminary test showed that PCs of all spectral bands can contribute to the estimation of canopy chlorophyll content and they were also used in previous studies for retrieving vegetation properties [69]. Therefore, top 10 PCs that can cover the majority of spectral band information (i.e., PC1 to PC10) were calculated and selected as predictor variables. Total 39 predictor variables were used to build LRs and their performances were later evaluated.

2.3.2. Partial Least Square Regression

The partial least square regression in Python with Scikit-Learn library was applied in this study [88]. The 29 VIs and 10 PCs used in LRs were also used as predictor variables in PLSR. In addition, since spectral bands may also contribute to the model prediction, reflectance values of the 301 bands were thus also used as predictor variables. These variables were named as Re + wavelength, such as Re550 and Re800. Moreover, image textural features (e.g., entropy, homogeneity) contain essential information describing vegetation structural patterns and spatial variations, which potentially can also contribute to the estimation of vegetation properties [55,89]. Therefore, eight textural variables, including mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment, and correlation, were calculated for selected spectral bands in ENVI using a kernel size of 5 × 5 (close to the size of a study site). These textural variables were named as b + wavelength–variable type, such as b800-mean and b740-entropy. Calculating eight textural variables for all of the spectral bands are computationally challenging with much information overlap since some spectral bands are correlated. Therefore, seven bands that are most frequently used for calculating the VIs listed in Table 1, were selected for calculating image textural variables. Details of these seven bands can be found in Section 3.1. Overall, there were 29 VIs, 10 PCs, 301 reflectance, and 56 textural variables used as predictor variables in the PLSR.

The PLSR is capable of dealing with a large number of predictor variables by converting them into a few latent variables. The predicted residual error sum of squares (PRESS) statistics was applied to identify the optimal number of latent variables in PLSR that has the lowest model error. Since some predictor variables may not contribute to the regression model or may bring in noise, a backward feature elimination method was applied to eliminate the less promising variables and identify the model that achieves the highest accuracy [19,48]. The variable importance on projection (VIP) values were calculated in PLSR to evaluate effects of variables on the regression model [90]. The least important variables (i.e., with lowest VIP values) were progressively deleted and the model was consequently rebuilt. Thus, a series of models were established and then also validated. The model with the lowest RMSE was regarded as the best performing model [48].

2.3.3. Random Forest Regression

The random forest regression in Python with Scikit-Learn library was utilized in this study [88]. There are three essential parameters need to be determined in the RFR model: number of trees to grow (Ntree, default 100), number of predictor variables to test at each splitting node (Mtry, default all variables), and node size that is the minimum number of samples required to be at a leaf node (Nodesize, default 1) [34,44]. The Mtry and Nodesize values were kept as default as suggested in previous studies, while the Ntree was tested using 100, 200, 500, 1000, and 2000 to optimize the model performance (i.e., lowest RMSE) [33,34,41,46].

The predictor variables used in PLSR, including VIs, reflectance, PCs, and textural variables, were also used in RFR. Importance values of all predictor variables were calculated in RFR [44], and then the backward feature elimination method was applied to RFR to progressively remove the least important variables [19,33,34,41]. The model with the lowest RMSE was regarded as the optimal model.

2.3.4. A Modified PROSAIL

The PROSAIL model was downloaded from this link (http://teledetection.ipgp.jussieu.fr/prosail/, accessed in June 2017, version 5B) and was operated in Matlab [11]. This original PROSAIL was found to suffer when simulating the spectra of heterogeneous canopies [22,53], likely because the leaf-level PROSPECT model cannot sufficiently simulate spectra of senescent leaves [54]. Therefore, a modified PROSPECT that is capable of tackling this problem was developed and applied in this study to generate a modified PROSAIL [54]. Parameters in PROSAIL include leaf chlorophyll content, leaf water content, leaf area index, solar and view geometry, and soil reflectance [9,11,54]. Vegetation canopy reflectance was extracted from hyperspectral imagery. A lookup table approach (LUT) was applied as the inversion method to estimate vegetation properties from the canopy reflectance data [50,91]. A brief workflow of the inversion is shown in Figure 4. A LUT was defined using all possible values within a value range for each parameter in the PROSAIL model. This value range was determined with prior knowledge from the literature [9,11,54]. For instance, the leaf chlorophyll content has a range of 0~100 µg/cm² and the leaf water content (i.e., equivalent water thickness) has a range of 0~0.05 cm. Then, the model calculated all possible combinations of values from different parameters and simulated reflectance correspondingly. This method covers all types of leaves and canopy structures, thus generating a dataset including all possible reflectance. With the measured canopy reflectance from hyperspectral imagery, a global search was performed in the LUT to identify parameter combinations that yield the best fit between measured and forward-simulated vegetation spectra. At last, vegetation properties, i.e., optimal parameter combinations, including leaf chlorophyll content and LAI, were determined.

3. Results and Discussion

3.1. Analysis of Predictor Variables

The correlations between each predictor variable (e.g., VIs, reflectance, PCs, and textural variables) and canopy chlorophyll content were first examined and the top 30 strongest ones (or groups of neighboring reflectance) are illustrated in Figure 5. These top 30 strongest correlations have Pearson’s r values around 0.8. The corresponding top 30 predictor variables include VIs, PCs, reflectance, and textural variables, indicating these different types of variables can all potentially contribute to the estimation of canopy chlorophyll. For the June image that captured homogeneous green vegetation, the top 17 variables are all VIs that have Pearson’s r values over 0.8, such as MCARI1, ZM, and mSR3 (Figure 5a). The VIs were selected based on their good performance in previous studies for estimating canopy chlorophyll and thus their strong correlations were expected in this study. Several other types of predictor variables, such as PC2, b800-Mean, and Re801–850, have only slightly weaker correlations with canopy chlorophyll than the VIs.

To further understand which bands were used in the VIs listed in Table 1 and identify the most frequently used ones that may be important for estimating canopy chlorophyll, the frequencies of different bands used in VIs were calculated and plotted in Figure 6. The top four most often used bands are centered around 550, 670, 750, and 800 nm, which are in the high reflection range in green, the deep absorption portion in red, the increasing reflection part in red-edge, and the peak reflection range in NIR, respectively (Figure 6). The reflection strength in green and red is influenced by leaf chlorophyll, while the strength in NIR is related to leaf structure and LAI (i.e., canopy density) [28]. Bands in the blue ranges, such as those around 450 nm, also have high frequencies. In addition, the bands in the red-edge range, such as those between 700 and 750 nm (e.g., 704 nm, Figure 6), are also often used in VIs. It is reported that using red-edge bands in VIs can minimize the saturation problem of the indices, reduce the influences of atmospheric and water absorption and background noise, and mitigate the vegetation surface scattering and bidirectional reflectance distribution function (BRDF) effects [33]. Several VIs that have used red-edge bands, such as MCARI1, ZM, and mSR3, show very strong correlations with the canopy chlorophyll (Figure 5). Overall, seven bands, including 450,550,670,680,704,750, and 800 nm, are more frequently used in the VIs, indicating their importance for estimating chlorophyll. Images of these seven bands are utilized for producing textural variables mentioned in Section 2.3.2.

For the August image that is featured by senescent canopies (i.e., with mixed green and senescent vegetation and low canopy chlorophyll), the VIs have weaker correlations with canopy chlorophyll than the other types of predictor variables (Figure 5b). For instance, the top five variables (or groups of variables) are reflectance (e.g., Re901–950) and PC1. Only three VIs (i.e., SPVI, DVI, and TSAVI) are among the top 10. This indicates that the VIs have weaker correlations with chlorophyll of senescent canopies than that of green canopies. Many VIs are built based on low reflectance in the red and high reflectance in the NIR (as green canopy reflectance shown in Figure 6) to highlight the green vegetation information and chlorophyll signal. However, for senescent canopies, the reflectance is higher in red and lower in NIR (as senescent canopy reflectance shown in Figure 6). The effectiveness of VIs for highlighting the chlorophyll information of senescent canopies is lower and thus they have relatively weaker correlations compared to that of green canopies. In contrast, reflectance in the NIR (e.g., 800–1000 nm) has strong correlations with chlorophyll of senescent canopies (Figure 5b). This is probably because the reflectance strength in NIR is highly affected by LAI and the canopy chlorophyll content was calculated using LAI. Canopies with high LAI show high reflectance in NIR and also have high canopy chlorophyll content, which thus leads to the strong correlation between reflectance in NIR and the canopy chlorophyll content. The results suggest that it is important to evaluate the influence of LAI and leaf chlorophyll content on the variations of canopy chlorophyll content. We found that the correlations between LAI and canopy chlorophyll are very strong (e.g., Pearson’s r of 0.96 for the June survey and 0.97 for the August survey), while that between leaf chlorophyll and canopy chlorophyll are much weaker (e.g., Pearson’s r of 0.31 for the June survey and 0.28 for the August survey). These results demonstrated that LAI has much stronger influence on the variation in canopy chlorophyll content than the leaf chlorophyll content.

3.2. Optimization of PLSR and RFR

The importance of each predictor variable was evaluated in the PLSR and RFR, respectively, and a backward feature elimination method was used to remove less important variables from the model with the aim of identifying the best performing PLSR and RFR (i.e., with the lowest RMSE). The selected variables for the optimal PLSR and RFR of the June or August image (hereafter named PLSR-June, RFR-June, PLSR-August, and RFR-August) are listed in

Table 2. Selected variables in the optimal partial least square regression (PLSR) and random forest regression (RFR) for June and August images, respectively. Variables are grouped by types (e.g., vegetation indices (VIs), PCs). If a large number of variables of one type were selected, only the top five (e.g., with the highest importance in the model) are listed as examples.

	June Image		August Image
	PLSR	RFR	PLSR	RFR
Selected Variables	VIs: mSR3 VREI MTCI MCARI1 REPI (Total 27)   Reflectance: Re628-Re1000 (Total 167)   PCs: PC1, PC2, PC4, PC5   Textural: b550-Homogeneity b670-Second Moment b704-Entropy b750-Mean b800-Mean (Total 16)	VIs: NDRE RENDVI ZM RRVI MRENDVI (Total 14)   Reflectance: Re502   PCs: PC5   Textural: b680-Mean	VIs: SPVI DVI TSAVI MSAVI RDVI (Total 11)   Reflectance: Re714-Re1000 (Total 144)   PCs: PC1   Textural: b800-Meanb 750-Mean	VIs: TSAVI DVI NDRE MCARI1 MSAVI (Total 7)   Reflectance: Re802-Re1000 (Total 24)   PCs: PC3   Textural: b800-Mean b680-Correlation

. Different types of variables, including VIs, reflectance, PCs, and textural variables, are all among the selected variables of different models, indicating they can all contribute to model predictions.

Comparing the number of selected variables in PLSR and RFR, the PLSR needs to select more variables than the RFR does to achieve the optimal performance. For instance, the PLSR-June selected 27 VIs,167 reflectance, 4 PCs, and 16 textural variables to reach the highest accuracy, while the RFR-June selected 14 VIs, 1 reflectance, 1 PC, and 1 textural variable to achieve the best performance. Similar results can be found for the PLSR-August and RFR-August. This is because the PLSR and the RFR use different techniques to evaluate variable importance in the model and then select the variables based on their importance. Specifically, the PLSR considers the predictive power of variables when converting them to latent variables [36,90]. Therefore, the PLSR selected a large number of variables to gain the maximum predicting power, even when these variables may be correlated (e.g., the 167 reflectance in PLSR-June and the 144 reflectance in PLSR-August are intercorrelated, respectively). In contrast, the RFR selects the best performing variable for splitting nodes and growing trees and rates it as more important in the model [17,46]. Other variables that are correlated with this one will have limited contribution to the model and thus are rated as less important. Therefore, compared to PLSR, a smaller number of predictor variables are needed in the RFR model to achieve optimal model performance.

3.3. Forward Simulation Using PROSAIL

A modified PROSAIL model, which integrated a modified PROSPECT-5 model and SAIL model, was used for simulating vegetation spectra (e.g., leaf and canopy reflectance) and estimating vegetation properties (e.g., leaf chlorophyll content, LAI). The modified PROSPECT-5 model worked well for simulating reflectance of different types of leaves (e.g., green and senescent), achieving an average RMSE of 0.009 with a standard deviation of 0.004. Examples of leaf spectral simulations are shown in Figure 7a,b. The modified PROSAIL model also performed well for simulating reflectance of both green and mixed canopies, achieving an average RMSE of 0.029 with a standard deviation of 0.018. Examples of canopy spectral simulations are shown in Figure 7c,d. The effectiveness of PROSAIL for estimating vegetation properties are evaluated in the next sections.

3.4. Result Comparison of Different Methods

LR, PLSR, RFR, and PROSAIL were utilized for estimating canopy chlorophyll content from images acquired in June and August, respectively. The validation results, including R² and RMSE, are shown in Figure 8. 39 LRs were built using 29 VIs and 10 PCs, and their accuracy values are described in boxplots in Figure 8. Overall, the best performing LRs, along with PLSR, RFR, and PROSAIL achieved good accuracies, such as R² ~ 0.82 and RMSE ~ 15.0 µg/cm² for the June image and R² ~ 0.75 and RMSE ~ 17.0 µg/cm² for the August image. The estimation accuracy for the August image is slightly lower than that of the June image, which is possibly due to vegetation senescence and reduced canopy chlorophyll content. Different models also have different performances for estimating canopy chlorophyll. For the June image, half of LRs achieved R² > 0.6 and RMSE < 20.0 µg/cm² (Figure 8a,b). A few top performing LRs included the ones using MCARI1 (R² = 0.80, RMSE = 13.9 µg/cm²), ZM (R² = 0.75, RMSE = 15.7 µg/cm²), and mSR3 (R² = 0.73, RMSE = 16.3 µg/cm²). For the August image, effectiveness of LRs are obviously lower than that for the June image. Only half of LRs achieved R² values higher than 0.4 and RMSE lower than 26.0 µg/cm² (Figure 8c,d). A few top performing LRs included the ones using PC1 (R² = 0.74, RMSE = 17.7 µg/cm²), SPVI (R² = 0.71, RMSE = 18.7 µg/cm²), and DVI (R² = 0.71, RMSE = 18.8 µg/cm²). The less efficient performance of LRs for the August image is probably due to the relatively low canopy chlorophyll content in the late growing season and the predictor variables (e.g., VIs) are not as sensitive to the canopy chlorophyll as that in June. While for the PLSR, RFR, and PROSAIL, the accuracies only vary slightly for the June or August images, indicating the higher stability of these three models over the LRs.

The PLSR achieved higher accuracies than the LRs, which is expected as it utilized multiple predictor variables in the model (Figure 8). However, its accuracy values are only slightly higher than that of the top-performing LRs. For instance, for the June image, the PLSR reached an R² of 0.84 and an RMSE of 12.5 µg/cm², while the LR using MCARI1 achieved an R² of 0.80 and an RMSE of 13.9 µg/cm² (Figure 8a,b). This indicates that the LRs with appropriate VIs can potentially achieve performances close to that of the advanced regression models (e.g., PLSR). When considering LRs for practical applications, it is thus likely best to compare a few VIs that performed well in previous studies and select the one with highest accuracy for building LR.

Hyperspectral information is critical for the good performance of LRs and PLSR. For instance, the PLSR achieved an R² of 0.77 and RMSE of 16.6 µg/cm² for the August image (Figure 8c,d). In contrast, in a previous study that was conducted in the same study area using multispectral image and PLSR [55], the accuracy of canopy chlorophyll estimation (R² = 0.31, RMSE = 18.6 µg/cm²) is much lower than that achieved in this study using the hyperspectral imagery. This is expected because the predictor variables extracted from the hyperspectral image (e.g., VIs and PCs) can provide more information than that extracted from the multispectral image. Previous studies have also suggested narrowband indices are more stable and sensitive to vegetation properties while broadband indices are more likely impacted by environmental factors (e.g., view geometry, atmospheric influences) [28,35,74,92].

The RFR performed the best, acquiring an R² of 0.86 and RMSE of 12.1 µg/cm² for the June image and an R² of 0.81 and RMSE of 15.4 µg/cm² for the August image (Figure 8). This is similar to the result in previous studies that RFR had a higher accuracy than that of other models for estimating vegetation properties [17,19]. The good performance of RFR is probably owning to its advantages in using randomly selected subset training data at each node and in selecting the best-performing predictor variable for splitting the node (thus not being sensitive to noise), as well as its ability to handle the collinearity problem [44]. However, random forest is a type of ‘black box’ since the tree structures are not clear and not easily understandable [44,45]. It is thus difficult to find how the variables are applied in the trees, which is a drawback of RFR.

The PROSAIL also performed well, achieving an R² of 0.80 and RMSE of 21.4 µg/cm² for the June image and an R² of 0.75 and RMSE of 20.1 µg/cm² for the August image (Figure 8). Compared to the regression models, the PROSAIL performed similarly well as the best performing LR, although not as well as the PLSR and RFR (Figure 8). One essential advantage of PROSAIL is that it can be transferred to different sites in different seasons (e.g., June and August in this study), while regression models are mostly not transferable (i.e., require collecting new training data). In addition, PROSAIL can generate a set of vegetation biophysical and biochemical properties at the same time, such as leaf chlorophyll content, leaf water content, and LAI, which can provide a more solid understanding of vegetation growth condition and physiological status. While the regression models, including LR, PLSR, and RFR, can only retrieve one vegetation property (i.e., one dependent variable) at a time.

Scatter plots of measured canopy chlorophyll contents against estimated values from different models that were built for the June image are shown in Figure 9 as examples, aiming to further evaluate performance of different models (e.g., over- or under-estimation). For the top performing LRs (e.g., R²~0.75, RMSE~15 µg/cm²), such as the ones built with MCARI1 or mSR3 (Figure 9a,b), no clear over- or under-estimation is observed. For the LRs established with TSAVI or OSAVI1 (Figure 9c,d), which have medium accuracy values (e.g., R² ~ 0.62, RMSE ~ 19 µg/cm²), under-estimation occurs for the measured chlorophyll higher than 100 µg/cm². This is probably due to the saturation problem of TSAVI and OSAVI1. These two indices are built with NIR and red bands (Table 1) and thus may suffer from the spectral saturation when the canopy is dense. Involving the red-edge bands, such as the top performing MCARI1 and mSR3, can probably mitigate this problem. The LRs built with GI or PPR have poor accuracy values (e.g., R² < 0.35, RMSE > 25 µg/cm²) and show clear over-estimation for measured chlorophyll lower than 70 µg/cm² and under-estimation when over 70 µg/cm² (Figure 9e,f). These two indices are generated with only visible bands, such as GI with green and red and PPR with green and blue. They cannot capture spectral variations in the NIR range that is sensitive to the canopy density and thus had poor performance. Therefore, when selecting VIs for building LRs, it is essential to select a variety of VIs that use bands in different spectral ranges (e.g., visible, red-edge, and NIR) and evaluate their performance. The PLSR and RFR used a wide range of VIs and thus did not suffer from over- or under-estimation (Figure 9g,h). The PROSAIL is a physically-based model and uses full spectrum information, and thus also did not experience obvious over- or under-estimation (Figure 9i).

4. Conclusions

LR, PLSR, RFR, and PROSAIL models were compared in this research for estimating vegetation canopy chlorophyll content from bi-seasonal hyperspectral images in a heterogeneous grassland area. Overall, the best performing LRs, along with PLSR, RFR, and PROSAIL, performed well in the estimation of canopy chlorophyll content, achieving R² ~ 0.80 and RMSE ~ 16.0 µg/cm². Accuracies of these four models in the August image is lower than that of the June image, especially for LRs, indicating the influence of vegetation phenological changes (e.g., vegetation senescence, increased heterogeneity, reduced chlorophyll content) on the performance of models. The PLSR and RFR had better performance than LRs and PROSAIL, while RFR is the best performing one with an R² of 0.86 and an RMSE of 12.1 µg/cm² for the June image and an R² of 0.81 and an RMSE of 15.4 µg/cm² for the August image. For ease of operation in practical projects for estimating canopy chlorophyll, one can use LRs by testing a couple of VIs that performed well in previous studies and selecting the one that has the best performance. It is also critical to select different VIs that are built with bands in different spectral ranges (e.g., visible, red edge, and NIR), since canopy chlorophyll is sensitive to different spectral ranges when it is high or low (i.e., different seasons). To achieve a higher accuracy and a more stable model, we suggest using PLSR or RFR. PROSAIL is a bit more complicated than the regression models and requires a set of parameters. However, the model itself is transferable and it can retrieve a wide range of vegetation biophysical and biochemical properties at the same time. PROSAIL is thus suggested if the model needs to be applied on different sites over different seasons or for retrieving a couple of vegetation properties. One should keep in mind that PROSAIL is not directly applicable to all ecosystems with various conditions (e.g., forests or crops, homogeneous or heterogeneous canopies). Appropriate model evaluations with necessary model modifications or parameter adjustments are strongly recommended.

Different types of image features, including VIs, reflectance, PCs, and textural variables, were applied in the PLSR and RFR as predictor variables and it was found that they can all contribute to model predictions (e.g., all the four types of variables were selected in the optimal models). Both PLSR and RFR are capable of dealing with a large number of variables and additional techniques (e.g., backward feature removal) can be applied to remove the less important variables and improve model efficiency and accuracy. The importance values of different variables evaluated by PLSR and RFR can help to understand their contributions to the models. The PLSR and RFR use different approaches to evaluate variable importance and thus the important variables that they selected are also different. The PLSR requires more variables than the RFR to achieve the optimal performance. RFR is a type of ‘black box’ since the tree structures are not clear and not easily understandable. Further research on understanding RFR results is warranted. For the images acquired in different seasons, the top important variables ranked by the models are different. As a result, variables selected in the optimal models are also different. This indicates that one predictor variable may have very different contributions to the models if using images acquired in different seasons.