Integration of Landsat-8 Thermal and Visible-Short Wave Infrared Data for Improving Prediction Accuracy of Forest Leaf Area Index

Abstract

Leaf area index (LAI) has been investigated in multiple studies, either by means of visible/near-infrared and shortwave-infrared or thermal infrared remotely sensed data, with various degrees of accuracy. However, it is not yet known how the integration of visible/near and shortwave-infrared and thermal infrared data affect estimates of LAI. In this study, we examined the utility of Landsat-8 thermal infrared data together with its spectral data from the visible/near and shortwave-infrared region to quantify the LAI of a mixed temperate forest in Germany. A field campaign was carried out in August 2015, in the Bavarian Forest National Park, concurrent with the time of the Landsat-8 overpass, and a number of forest structural parameters, including LAI and proportion of vegetation cover, were measured for 37 plots. A normalised difference vegetation index threshold method was applied to calculate land surface emissivity and land surface temperature and their relations to LAI were investigated. Next, the relation between LAI and eight commonly used vegetation indices were examined using the visible/near-infrared and shortwave-infrared remote sensing data. Finally, the artificial neural network was used to predict the LAI using: (i) reflectance data from the Landsat-8 operational land imager (OLI) sensor; (ii) reflectance data from the OLI sensor and the land surface emissivity; and (iii) reflectance data from the OLI sensor and land surface temperature. A stronger relationship was observed between LAI and land surface emissivity compared to that between LAI and land surface temperature. In general, LAI was predicted with relatively low accuracy by means of the vegetation indices. Among the studied vegetation indices, the modified vegetation index had the highest accuracy for LAI prediction (R²_CV = 0.33, RMSE_CV = 1.21 m²m⁻²). Nevertheless, using the visible/near-infrared and shortwave-infrared spectral data in the artificial neural network, the prediction accuracy of LAI increased (R²_CV = 0.58, RMSE_CV = 0.83 m²m⁻²). The integration of reflectance and land surface emissivity significantly improved the prediction accuracy of the LAI (R²_CV = 0.81, RMSE_CV = 0.63 m²m⁻²). For the first time, our results demonstrate that the combination of Landsat-8 reflectance spectral data from the visible/near-infrared and shortwave-infrared domain and thermal infrared data can boost the estimation accuracy of the LAI in a forest ecosystem. This finding has implication for the prediction of other vegetation biophysical, or possibly biochemical variables using thermal infrared satellite remote sensing data, as well as regional mapping of LAI when coupled with a canopy radiative transfer model.

1. Introduction

Leaf area index (LAI) is extensively applied to observe and monitor ecosystem functions (e.g., vegetation growth, and physiological activity) [1,2,3]. Due to the control of LAI over primary production (e.g., photosynthesis), transpiration, evapotranspiration, energy exchange as well as other physiological characteristics pertinent to the wide range of ecosystem processes, the accurate prediction of LAI has been a concern for a broad spectrum of studies [4,5,6,7,8,9]. Moreover, LAI has lately been suggested as being one of the essential biodiversity variables (EBVs) that are suitable for satellite monitoring, among many other variables [10,11]. In addition, the demand for LAI monitoring over large areas in recent years has increased due to the monitoring and modelling of climate change as well as habitat degradation [12,13].

LAI has been widely retrieved by means of visible/near-infrared (VNIR, 0.3–1.0 µm) as well as shortwave-infrared (SWIR, 1.0–2.5 µm) spectral data over different ecosystems with varying degree of success [14,15,16,17,18,19,20,21]. In this respect, Verrelst et al. [22] were evaluated the all possible band combinations for two- and tree-bands indices as well as different machine learning approaches using Sentinel-2 data for LAI retrieval and revealed machine learning approaches performed with greater accuracy. Additionally, Neinavaz et al. [23] demonstrated that multivariate methods (e.g., artificial neural network) are the most promising approach in comparison with the univariate approaches (e.g., vegetation indices) for prediction of the LAI using hyperspectral thermal infrared (TIR, 8–14 µm) data. However, the potential of TIR remote sensing for estimating vegetation biophysical variables in general, and LAI in particular, has not been sufficiently studied.

The advantage of TIR data in remote sensing vegetation studies is that homonuclear diatomic molecules (e.g., N₂ or O₂) do not demonstrate substantial spectral features in wavelengths between 8 and 14 µm; consequently, the atmosphere is mostly transparent over this domain [24]. The considerable advantage of using TIR data regarding LAI prediction despite the importance, maturity, and availability of multispectral remotely sensed data in the VNIR/SWIR domain is that saturation does not occur, not even at comparatively high LAI values [25]. Recently, under controlled laboratory conditions, LAI has been successfully quantified by means of TIR emissivity spectra [23]. However, the applicability and transferability of laboratory studies to LAI estimation, at the landscape level, using TIR data (i.e., land surface temperature (LST) and land surface emissivity (LSE) [26]), remain to be elucidated. LSE is a measure of the inherent energy of the surface in transforming kinetic energy into radiant energy over the surface [27]. Knowing the LSE is necessary for estimating the energy budget, evapotranspiration, water, and energy balance [28,29,30].

In addition, the LSE is a critical parameter for retrieving the LST [31]. When only limited thermal data (i.e., one thermal band) is available, LSE becomes even more critical for estimating LST [32,33]. LST is also an important parameter for describing the physical processes of surface energy and the hydrological cycle, and serve to broaden our knowledge on the surface equilibrium state concerning temporal and spatial changes [34]. As saturation has been observed in the prediction of LAI over the VNIR/SWIR domain using multispectral data, particularly in the forest ecosystem [3,35,36,37], as well as hyperspectral data at high LAI values [14], the integration of the VNIR/SWIR and TIR data may address this issue and boost LAI retrieval accuracy, as saturation does not occur at relatively high LAI values by means of TIR hyperspectral data [23]. In this regard, Breunig et al. [38] showed that a combination of SWIR (i.e., reflectance) and TIR (i.e., emissivity) data enhanced the discrimination between exposed soil and non-photosynthetic vegetation. More recently, Mushore et al. [39] found that the integration of Landsat-8 data from the thermal infrared sensor (TIRS, 10.6–12.5 µm) and the operational land imager (OLI, 0.4–2.5 μm) sensor achieved significantly higher accuracy when classifying urban landscapes. More recently, Bayat et al. [40] demonstrated that the satellite optical and TIR data could be successfully integrated to capture drought effects at the canopy level.

To the best of our knowledge, the potential for utilising a combination of Landsat-8 VNIR/SWIR and TIR data (e.g., LST and LSE) has not been investigated in the context of retrieving forest biophysical parameters, such as LAI. As little research has been done in this direction, we addressed this gap in our study. The aim of this study was dedicated to improving forest LAI retrieval accuracy by integrating Landsat-8 OLI and TIRS data through the calculation of LSE, and LST.

2. Material and Methods

2.1. General Description of the Study Area

Field measurements were performed over the Bavarian Forest National Park (BFNP), which is located in the federal state of Bayern, in the southeastern part of Germany, along the border with the Czech Republic (49°3′19″N, 13°12′9″E) (Figure 1). In total, the BFNP area is 24,250 hectares [41]. The elevation in this area ranges from 600 to 1453 m above the sea level. The BFNP has a temperate climate, and precipitation varies from 1200 to 1800 mm/year (of which 50% is snow). In some years, the highest rainfall is more than 2000 mm [41]. The minimum average annual temperature is between 3 °C and 6 °C. Three major forest types are recognisable in the BFNP. At the highest elevations (i.e., above 1100 m), Norway spruce (Picea abies) occurs with sub-alpine Spruce and a few Mountain ash (Sorbus aucuparia); At altitudes that encompass the slopes (i.e., between 600 and 1100 m), there are Norway spruce, Silver fir (Abies alba), European beech (Fagus sylvatica) and Norway maple (Acer pseudoplatanus). In wet depressions in the valleys, spruce forests mingled with mountain ash, Norway spruce, as well as Birches (Betula pendula, and B. pubescens) can be found [42,43,44]. In general, the dominant tree species in the BFNP are Norway spruce (67%), European beech (24.5%), and fir (2.6%) [45].

2.2. Collection of In Situ Structural Canopy Parameters

A field campaign was conducted in August 2015. The BFNP is covered in broadleaf, needle leaf (conifer) as well as mixed forest stands. Random sampling was chosen as the most straightforward strategy, and 37 plots were selected (plot size: 30 m × 30 m), resulting in four broadleaves, 26 needle leaves, and seven mixed forest plots. To record the coordinates of the centre of each plot, a Leica GPS 1200 system (Leica Geosystems AG, Heerbrugg, Switzerland) was used, achieving a roughly 1 m accuracy after post-processing [46]. For each plot, the plants species were determined, and the proportion of vegetation cover (P_V) and LAI, representing the structural forest parameters, were computed. LAI is a dimensionless parameter which is defined as the total of the one-sided leaf area (m²) per unit horizontal surface area (m²) [47]. A plant canopy analyser (LAI-2200, LICOR Inc., Lincoln, NE, USA) was used for measuring LAI in the field. Reference samples of the above-canopy radiation for each plot were collected through quantifying the incoming radiation in nearby open spaces under cloud-free conditions. Eventually, five below-canopy samples were quantified, and the LAI value was then computed by averaging these measurements for each plot. P_V has also been termed as the fractional vegetation cover and was initially introduced by Deardorff [48] as the proportion of the vertical projection area of the plant on the surface of the ground (including leaves, stalks, and branches) to the total vegetation area [49]. Using this definition, the P_V of each plot was computed using five upward-pointing digital hemispherical photographs (DHP), following Zhou et al. [50]. For each plot, five upward-pointing DHPs were collected. The images were acquired using a Canon EOS 5D, equipped with a fish-eye lens (Sigma 8 mm F3.5 EX DF), levelled on a tripod at approximately breast height (1.3 m above the ground) [51]. Each image had a high resolution of 5600 × 3898 pixels. The two-corner classification procedure [52] was used to minimise subjective thresholding on the blue channel of all obtained images so as to classify sky and canopy pixels. Further, the CanEye software was applied to estimate P_V by importing binary classified images. The arithmetic mean of P_V estimated from the five images was then considered to be the P_V of each plot. P_V was also used to calculate the LSE.

2.3. Satellite Data and Processing

The Landsat-8 data were acquired on 9 August 2015 for the study area (Landsat-8 Scene ID = LC81920262015221LGN01). The Landsat-8 satellite has two main sensors, the OLI and the TIRS (

Table 1. The Landsat-8 sensors, the operational land imager (OLI) and the thermal infrared sensor (TIRS) spectral bands, and their spatial resolution.

Landsat-8 Sensor	Bands	Wavelength (µm)	Resolution (m)
OLI	Band 1	0.43–0.45	30
	Band 2	0.45–0.51	30
	Band 3	0.53–0.59	30
	Band 4	0.64–0.67	30
	Band 5	0.85–0.88	30
	Band 6	1.57–1.65	30
	Band 7	2.11–2.29	30
TIRS	Band 10	10.60–11.19	100
	Band 11	11.50–12.51	100

). Since the Landsat-8 level-1 products were not atmospherically corrected, the OLI and TIRS images were corrected by converting digital numbers to radiance values, using coefficients supplied by the United States Geological Survey (USGS, https://landsat.usgs.gov/using-usgs-landsat-8-product). For the OLI images, the conversions of radiance to reflectance and atmospheric correction have been done using the FLAASH module, as the FLAASH properties consider water vapor, distribution of aerosols as well as scene visibility for the atmospheric correction. In this study, the normalised difference vegetation index threshold method (NDVI^THM) was applied to estimate LSE and LST. Therefore, the atmospheric correction was not needed for the TIRS bands [53]. At the outset, TIRS bands were acquired at 100 m resolution, but since 2010, they have been resampled to 30 m by the USGS, using cubic convolution to match with the OLI spectral bands. As such, the parallel use of Landsat-8 spectral and thermal data has been investigated in several studies [54,55,56].

2.4. Land Surface Emissivity and Land Surface Temperature

Several approaches exist for estimating LSE and LST using remotely sensed data [57]. Among these approaches, the NDVI^THM [58], which has been further modified and developed by Sobrino and Raissouni [59] and Valor and Caselles [60], has been considered to be a practical approach [27]. In the NDVI^THM approach, the statistical relationship between NDVI and emissivity over TIR spectral bands is used to determine LSE. In this respect, we derived LSE and LST using Landsat-8 images to retrieve the LAI for the BFNP. Band 10 of Landsat-8 was considered because instability in the calibration of band 11 has been reported by Barsi et al. [61]. The LSE can be computed though the relationship between the NDVI and the vegetation and soil emissivity as follows [57,59]:

$$LSE=\{\begin{array}{}\mathrm{(1a)}& NDVI<0.2& a_{\lambda }+b_{\lambda }{\rho }_{red}\mathrm{(1b)}& NDVI\ge 0.5& {\mathsf{\epsilon }}_{v\lambda }+d_{\mathsf{\epsilon }}\mathrm{(1c)}& 0.2\le NDVI\le 0.5& {\mathsf{\epsilon }}_{V\lambda }{P}_V+{\mathsf{\epsilon }}_{S\lambda }\times \left(1-P_V\right)+d_{\mathsf{\epsilon }}\end{array}$$

where $a_{\mathsf{\lambda }}$ and $b_{\mathsf{\lambda }}$ are channel-dependent regression coefficients, ${\rho }_{red}$ is a reflectivity value in the red region, and ${\mathsf{\epsilon }}_{v\lambda }$ and ${\mathsf{\epsilon }}_{s\lambda }$ are TIR band emissivity values for vegetation and bare soil, respectively. Both ${\mathsf{\epsilon }}_{v\lambda }$ and ${\mathsf{\epsilon }}_{s\lambda }$ can be measured directly in the field or downloaded from emissivity spectral libraries or databases. In this study, ${\mathsf{\epsilon }}_{v\lambda }$ and ${\mathsf{\epsilon }}_{s\lambda }$ were extracted from the MODIS University of California- Santa Barbara (USA) [62]. While P_V denotes the proportion of vegetation cover, dε stands for the cavity effect. Regarding flat surfaces, dε is inconsequential; however, for diverse and rough surfaces such as a forest ecosystem, dε can gain a value of 2% [63,64]. In addition, dε can be calculated by applying the following equation:

$$d_{\mathsf{\epsilon }}=\left(1-{\mathsf{\epsilon }}_s\right)\left(1-{{\rm P}}_v\right)F{\mathsf{\epsilon }}_v$$

(2)

where F is a shape factor, the mean value of which, assuming different geometrical distributions, is 0.55 [33,63].

To derive LST using TIR remotely sensed data, the brightness temperature (BT) should first be calculated, by means of the spectral radiance of TIRS bands, using the thermal constants [65]:

$$BT=\frac{K_2}{In \left(\frac{K_1}{L_{\lambda }}+1\right)}$$

(3)

where $L_{\lambda }$ denotes spectral radiance at the top of the atmosphere, and K₁ and K₂ are bands-specific thermal conversion constants, which are available from the metadata file of the Landsat-8 image. LST was calculated using the following equation, proposed by Stathopoulou and Cartalis [66]:

$$LST=BT/1+W\times \left(BT/p\right)\times Ln \left({\mathsf{\epsilon }}_i\right)$$

(4)

where W is the wavelength of emitted radiance, ${\mathsf{\epsilon }}_i$ denotes LSE and $p$ is equal to 1.438 × 10⁻² mK which is calculated with the following formula:

$$P=h\times C/S$$

(5)

where $h$ stands for the Planck’s constant (6.626 × 10⁻³⁴ Js), $S$ is the Boltzmann constant (1.38 × 10⁻²³ J/K) and $C$ denotes the speed of light (2.998 × 108 m/s). The LST is converted to Celsius by subtracting from 273.15, since it was derived in Kelvin.

2.5. Estimation of Leaf Area Index

2.5.1. Estimation of Leaf Area Index Using Vegetation Indices

Pearson’s r correlation was applied to determine the correlation between the LSE and LST and the in situ measured LAI. The most popular vegetation indices for retrieving vegetation parameters from VNIR/SWIR remotely sensed data are ratio-based [67]. In this study, eight different vegetation indices, which have been widely applied in the literature to derive the LAI, were examined (

Table 2. Vegetation indices which have been considered in this research concerning prediction of leaf area index.

Spectral Index	Original Equation	Abbreviation	Reference
Ratio Vegetation Index	$\ast \frac{{\rho }_{NIR}}{{\rho }_{Red}}$	SR	[69]
Modified Simple Ratio	$\frac{{\rho }_{NIR}-{\rho }_{Red}-1}{{\left({\rho }_{NIR}/{\rho }_{Red}\right)}^{0.5}+1}$	MSR	[70]
Difference Vegetation Index	${\rho }_{Red}-{\rho }_{NIR}$	SD	[71]
Renormalized Difference Index	$\frac{{\rho }_{NIR}-{\rho }_{Red}}{\sqrt{{\rho }_{NIR}+{\rho }_{Red}}}$	RDI	[72]
Modified Vegetation Index	$\frac{{\rho }_{NIR}-1.2{\rho }_{Red}}{{\rho }_{NIR}+{\rho }_{Red}}$	MVI	[73]
Normalized Difference Vegetation Index	$\frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+{\rho }_{Red}}$	NDVI	[74]
Enhanced Vegetation Index	$\ast \ast G\frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+C_1\times {\rho }_{Red}-C_2\times {\rho }_{Blue}+L}$	EVI	[75]
Reduced Simple Ratio	$\ast \ast \ast \frac{{\rho }_{NIR}}{{\rho }_{Red}}\left(1-\frac{{\rho }_{SWIR}-{\rho }_{SWI{R}_{min}}}{{\rho }_{SWI{R}_{max}}-{\rho }_{SWI{R}_{min}}}\right)$	RSR	[18]

* Where $\rho$ denotes the reflectance value at given wavelength λ, NIR represents the near-infrared reflectance; ** where G is the gain factor, C₁ and C₂ stand as the coefficients of the aerosol resistance term and L denotes the soil adjustment factor; *** ${\rho }_{SWIR}$ represents the near-infrared reflectance, and ${\rho }_{SWI{R}_{min}}$ and ${\rho }_{SWI{R}_{maz}}$ are the minimum and maximum SWIR reflectance values found in each image, respectively.

). The coefficient of determination (R²) between each index and LAI was considered to assess the strength of the relationship between the LAI and the proposed indices. The reliability of the model in estimating LAI was evaluated using a cross-validation procedure [68]. The cross-validated coefficient of determination (R²_CV) and cross-validated root mean squared error (RMSE_CV) were applied to assess the estimated LAI. All analyses were computed using MATLAB R2017b (Mathwork, Inc.).

2.5.2. Estimation of Leaf Area Index Using Artificial Neural Networks

The artificial neural network ANN consists of different layers including inputs, hidden layers as well as outputs. In this study, three scenarios were considered as input layers to the artificial neural networks [76] for LAI estimation. These included the reflectances of bands 1 to 7 from the OLI sensor (VNIR/SWIR), the combination of reflectance and LSE as well as the combination of reflectance and LST (

Table 3. Different inputs for estimation of leaf area index using the artificial neural network. LST and LSE represent the land surface temperature and land surface emissivity, respectively. OLI bands represent reflectance measurements from seven bands over visible–near-infrared and shortwave-infrared regions.

OLI Bands	TIRS Band	Input	Output
✓	-	7 bands	1
✓	✓ (i.e., LST)	8 bands	1
✓	✓ (i.e., LSE)	8 bands	1

). For network training, the Levenberg–Marquardt algorithm was used as the common training algorithm in backpropagation networks to develop models for LAI prediction. Atkinson and Tatnall [77] suggested that by raising the number of hidden layers, the network could tackle more complex datasets. However, still, there is no specific rule for defining the optimal number of hidden layer. Since the prediction accuracy of ANN is related to the number of neurons in the hidden layer, the ideal ANN size was determined by examining various numbers of the neurons. In this respect, the early stopping approach was used to avoid over-fitting. In this approach, the training network pauses as soon as performance on the validation dataset begins to worsen [78]. Linear regression analyses were performed between the retrieved and measured LAIs. As in Section 2.5.1, the reliability of the ANNs in retrieving LAI was evaluated using a cross-validation procedure and the predicted LAI was determined using the R²_CV and RMSE_CV. All calculations were performed in MATLAB R2017b (Mathwork, Inc.).

3. Results

3.1. Leaf Area Index and Proportion of Vegetation Cover

Summary statistics of the in situ measurements of LAI and P_V are presented in

Table 4. Summary statistics of the in situ measured proportion of vegetation cover (P_V) and leaf area index (LAI) over the Bavarian Forest National Park (n = 37).

Variables	Minimum	Maximum	Mean		Std. Deviation	Coefficient of Variation
			Statistic	Std. Error
LAI	0.50	5.86	3.34	0.24	1.46	43.60
PV	0.39	0.82	0.61	0.02	0.12	20.51

. The measured LAI from 37 sample plots demonstrated a maximum value of 5.86 m²m⁻². The measured P_V ranged from 0.39 to 0.82 (range = 0.43) and showed a maximum value of 0.82. The range of LAI values differed according to the forest stand. For instance, the LAI value of plots for broadleaf (mean = 3.04) and needle leaf (mean = 3.67) ranged from 0.51 to 4.89 and from 2.29 to 5.18, respectively, while the LAI values for the mixed forest stand was from 4.25 to 5.86 (mean = 4.77).

3.2. Relationships among Leaf Area Index, Land Surface Temperature, and Land Surface Emissivity

Applying the procedure defined for the estimation of LSE and LST in Section 2.4 to the Landsat-8 data, first, the LSE and LST were calculated, and then, their values for each plot were extracted over the BFNP. The relations between the LSE and LST obtained for each plot, and the corresponding LAI was then studied. A Pearson correlation coefficient revealed significant correlation between LSE and LAI (r = 0.66; P-value < 0.001). However, there was no significant correlation observed between LST and LAI (r = 0.26; P-value = 0.11). The relationships between LAI and LSE as well as LST for different plots, using a first-order polynomial, are presented in Figure 2. For the LSE, the regression parameters were negligible, with a slope of 0.00 and an offset of 0.98, while for the LST, the slope and offset were 1.13 and 22.21, respectively.

3.3. Estimated Leaf Area Index Using Vegetation Indices

LAI was predicted with moderate accuracy using the considered vegetation indices. Comparison of the R²_CV and RMSE_CV values among the investigated indices revealed that MVI retrieved an LAI with slightly greater accuracy compared with other vegetation indices (

Table 5. The coefficients of determination (R²) and cross-validation procedure among different indices calculated using reflectance over optical domain and leaf area index.

Vegetation Index	R2	Cross-Validation Procedure
		R2CV	RMSECV
SD	0.165	0.100	1.413
SR	0.373	0.308	1.230
RDI	0.234	0.166	1.357
MSR	0.292	0.227	1.305
MVI	0.408	0.331	1.218
NDVI	0.313	0.321	1.288
EVI	0.216	0.210	1.313
RSR	0.209	0.126	1.392

), and a linear relationship existed between the estimated and measured LAIs.

3.4. Estimating Leaf Area Index Using Artificial Neural Networks

Next, LAI was retrieved with a combination of spectral data from the VNIR/SWIR (i.e., reflectance) and the TIR data (i.e., LSE and LST, separately) using the ANN approach. Comparing the three scenarios that are presented in Table 3, the combination of spectral information from reflectance and LSE improved the prediction accuracy yielding an RMSE_CV of 0.63 m²m⁻², compared to using only the OLI data (RMSE_CV = 0.83 m²m⁻²) and a combination of the reflectance and LST (RMSE_CV = 0.63 m²m⁻²). The relationships between estimated and measured LAIs using the ANN model for different scenarios are shown in Figure 3. As can be seen, the prediction accuracy was increased in comparison with the vegetation indices. However, our findings suggest that the ANN calculated from the reflectance tended to overestimate the LAI values, which were less than 1 m²m⁻², whereas they performed with higher accuracy for LAI values between 2 m²m⁻² and 5 m²m⁻².

4. Discussion

For the first time, this study attempted to highlight the importance of integration of the satellite TIR and VNIR/SWIR data for improving the estimation of LAI over the temperate forest. The results of this study demonstrate that integration of the VNIR/SWIR and TIR satellite data has a high potential for boosting the retrieval accuracy of the LAI as the most important vegetation biophysical variable as well as the EBV. However, the LAI could be predicted with higher accuracy by using an integration of reflectance and LSE, calculated from TIRS band 10, rather than LST and reflectance. This observation can probably be explained by the fact that LSE is sensitive to LAI variations [25], and also is considered to be an indicator of material composition [27,53], while LST is a function of soil water content, surface soil, and the percentage of an area covered by vegetation [79], and so might be influenced by environmental conditions.

A review of the literature has shown that hyperspectral data are more efficient at providing extra information in comparison with multispectral data in quantifying vegetation characteristics over the VNIR/SWIR and TIR domains [80,81,82,83]. As in broadband sensors, the available information is usually masked [84,85]. Therefore, the prediction accuracy of LAI could be increased by applying VNIR/SWIR and TIR hyperspectral data.

An early paper by Neinavaz, Skidmore, Darvishzadeh and Groen [23] revealed that the highest accuracy in predicting LAI was obtained using the 10–12 µm wavelength range in combination with the bands in the 8–11 µm range. Therefore, having only one TIR band (i.e., the band 10 TIRS sensor) in the atmospheric window between 10 and 12 μm for Landsat-8 may have reduced the prediction accuracy of LAI. Moreover, the combination of several bands from different regions demonstrates better sensitivity than using one or two bands from a particular part of the electromagnetic spectrum for predicting vegetation parameters, since each band contains relevant and efficient information regarding its regions [86,87].

In this study, LAI was estimated with relatively moderate accuracy using vegetation indices (Table 5). However, the results further confirmed that LAI is predicted with higher accuracy using the ANN approach (Figure 3) than with vegetation indices. This result is in agreement with the findings of Danson et al. [88] and Neinavaz, Skidmore, Darvishzadeh and Groen [23], which respectively showed that LAI was successfully predicted using an ANN approach when either VNIR/SWIR or TIR hyperspectral data were utilised. The results of this study revealed that a combination of both TIRS (i.e., LSE) and reflectance data using a trained ANN (R²_CV = 0.81, RMSE_CV = 0.63 m²m⁻²), is more reliable, and achieved higher prediction accuracy than the use of reflectance and its combination with LST. It should be highlighted that in this study, the in situ P_V measurements were used for computing LSE. Hence, an accurate estimation of the P_V by means of remote sensing data for calculating LSE using NDVI^THM should be taken into account in future studies as a review of the literature has shown that P_V could be estimated with different accuracy degrees by means of vegetation indices [89] and machine learning approaches [90] over different ecosystems.

This study, besides showing that LAI is predictable using a combination of TIR and VNIR/SWIR data, also showed that ANN approaches improve model accuracy compared to univariate techniques for estimating vegetation characteristics, and have significant potential for the operational retrieval of LAI from remote sensing data [14,23,88].

5. Conclusions

The potential to predict LAI using integrated TIR and VNIR/SWIR data in a mixed temperate forest such as the BFNP has previously not been investigated. In this study, we have demonstrated this capability with the use of an ANN approach. The results demonstrate that the integration of LSE and reflectance data can improve LAI estimation. However, while the ANN-based models serve as a reliable approach to estimating LAI, this needs to be explored using hyperspectral satellite or airborne data over the TIR region. A linear relationship was found between the predicted and measured LAIs in all considered models. Our findings suggest that the combination of the VNIR/SWIR and TIR (i.e., LSE) data improved the prediction accuracy for estimating LAI, even under non-controlled conditions, which is a novel discovery. This once again reveals that TIR remotely sensed data are a valuable application in vegetation remote sensing studies. Further research is also required in the use of TIR hyperspectral data (i.e., airborne or spaceborne data), coupled with a canopy radiative transfer model. Our study offers practical techniques for estimating LAI as an important EBV, which will be valuable to the assessment as well as monitoring of biodiversity and ecosystem services.