Improving Field-Scale Wheat LAI Retrieval Based on UAV Remote-Sensing Observations and Optimized VI-LUTs

Abstract

Leaf area index (LAI) is a key biophysical parameter for monitoring crop growth status, predicting crop yield, and quantifying crop variability in agronomic applications. Mapping the LAI at the field scale using multispectral cameras onboard unmanned aerial vehicles (UAVs) is a promising precision-agriculture application with specific requirements: The LAI retrieval method should be (1) robust so that crop LAI can be estimated with similar accuracy and (2) easy to use so that it can be applied to the adjustment of field management practices. In this study, three UAV remote-sensing missions (UAVs with Micasense RedEdge-M and Cubert S185 cameras) were carried out over six experimental plots from 2018 to 2019 to investigate the performance of reflectance-based lookup tables (LUTs) and vegetation index (VI)-based LUTs generated from the PROSAIL model for wheat LAI retrieval. The effects of the central wavelengths and bandwidths for the VI calculations on the LAI retrieval were further examined. We found that the VI-LUT strategy was more robust and accurate than the reflectance-LUT strategy. The differences in the LAI retrieval accuracy among the four VI-LUTs were small, although the improved modified chlorophyll absorption ratio index-lookup table (MCARI2-LUT) and normalized difference vegetation index-lookup table (NDVI-LUT) performed slightly better. We also found that both of the central wavelengths and bandwidths of the VIs had effects on the LAI retrieval. The VI-LUTs with optimized central wavelengths (red = 612 nm, near-infrared (NIR) = 756 nm) and narrow bandwidths (~4 nm) improved the wheat LAI retrieval accuracy (R² ≥ 0.75). The results of this study provide an alternative method for retrieving crop LAI, which is robust and easy use for precision-agriculture applications and may be helpful for designing UAV multispectral cameras for agricultural monitoring.

1. Introduction

Crop growth varies between/within fields owing to differences in environmental conditions (e.g., climate, soil, and water), management (e.g., fertilization, irrigation, and pest/disease control) [1], and crop genotypes. Mapping the spatial variation of crop canopy traits and the factors affecting them is important for agriculture since these field-scale maps are helpful for precise agronomic management. The precise management is essential to improve crop production and light-utilization efficiency [2]. Satellite remote sensing has been widely applied in crop monitoring. However, its coarse spatial resolutions (e.g., moderate-resolution imaging spectroradiometer (MODIS) (1 km), Sentinel 1A (10 m), Landsat (30 m)) fail to fulfill the requirement of precise mapping at the field scale (~cm level) [3]. In contrast to satellite platforms, ground platforms are able to map at the centimeter scale, but they fail to be applied in relatively large areas [4]. However, unmanned aerial vehicles (UAVs) are able to solve these issues. UAVs are easy to use, have a low cost, and are able to capture data with a high spatial resolution (~cm); they provide a promising platform for monitoring the crop growth status and evaluating the results of field management. Multispectral cameras onboard UAVs capture the spectral reflectance of the vegetation canopy at specific bands. These remote-sensing observations can be used to obtain the leaf area index (LAI) which is a crucial biophysical parameter related to crop growth [3,5]. LAI data can be used to quantify crop variability [6], monitor crop growth [7], predict the crop yield [8], estimate the amount of aboveground biomass [9], and evaluate the effects of field management [10]. However, LAI retrieval methods based on UAV remotely sensed data must be robust and easy use so that the crop LAI will be estimated with similar accuracy to fulfill the requirement of practical application.

The PROSAIL, a radiative transfer model, which is run in a forward mode with appropriate parameter combinations (e.g., LAI, chlorophyll content, and leaf structure parameters) to simulate the spectral reflectance of vegetation canopy, is widely employed to retrieve LAI. The PROSAIL inputs (parameter combinations) and outputs (spectral reflectance) are used to generate look-up tables (LUTs) [11]. Thus, LAI, one of the PROSAIL inputs, can be retrieved by searching for the corresponding spectral reflectance in a given LUT using cost functions. The cost function is developed by simulated reflectance and measure reflectance. This LUT inversion method can reduce computational loads, lower the risk of converging to local minima, and help to avoid the ill-posed problems caused by uncertainties of the measurements and models [12].

The reflectance-LUT strategy, searching for the optimal simulated LAI by simulated reflectance and measured reflectance, is widely employed, particularly for traditional satellite remote sensing (Table 1). For UAV/airborne remote sensing, reflectance-LUTs are generated via multispectral or hyperspectral imaging. LAI retrieval based on UAV/airborne hyperspectral reflectance-LUTs can achieve high accuracy [13]. However, the multispectral cameras onboard UAVs typically have only three to six channels [14]. Spectral information for crops captured by multispectral cameras may be deficient owing to the limited bands. Therefore, the direct use of multispectral reflectance for LAI retrieval may be inaccurate. Generally, for multispectral UAV/airborne data, LAI estimation requires auxiliary methods, such as multi-angle observations for more vegetation information [12] and combinations with other algorithms (e.g., Bayesian networks [15], artificial neural networks, and random forests [16]). Multi-angle observations cost extra time and labor; machine learning requires ground sampling for the development of models. Therefore, these auxiliary methods may increase the difficulty of UAV multispectral remote sensing in practical application.

Vegetation indices (VIs) are calculated by different band combinations and mathematical formulas. Compared with spectral reflectance, vegetation indices (VIs) can provide more vegetation information, since VIs can offset parts of the errors from the surrounding environment and reveal implied plant information. Therefore, the VI-LUT strategy, searching for the optimal simulated LAI in a given LUT by simulated VIs (calculated from the PROSAIL outputs) and measured VIs, is an alternative strategy. Most of the previous studies related to VI-LUTs were based on satellite remote sensing [17]; no references were based on UAV/airborne remote sensing. Thus, we attempted to bridge this gap by investigating the robustness of the VI-LUT strategy for LAI retrieval using UAV multispectral remote sensing.

Table 1. Summary of LAI retrieval based on LUTs generated from the PROSAIL model.

LUTs	Platforms	Sensors	References
Reflectance-LUTs	UAV/Airborne	Multispectral camera	[12,15]
		Hyperspectral camera	[11,13]
	Satellite	ZY-3 MUX, GF-1 WFV, HJ-1 CCD	[18]
		Landsat	[19,20,21]
		Sentinel-2	[22,23,24]
		RapidEye	[22]
		CHRIS/PROBA	[25]
		MODIS	[20,26]
		AWiFS	[27]
	Ground-based	Hyperspectral spectrometer	[28,29,30]
VI-LUTs	satellite	CHRIS/PROBA	[31]
		RapidEye	[17]
		Landsat	[32]

LAI = Leaf area index; UAV = unmanned aerial vehicle; Reflectance-LUTs = reflectance-based lookup tables; VI-LUTs = Vegetation index-based lookup tables; ZY-3 MUX = the resources satellite No. 3 (ZY-3) with a multi-spectral camera (MUX); GF-1 WFV = the satellite GaoFen No. 1 (GF-1) with a wide field of view (WFV) camera; HJ-1 CCD = the environment satellite (HJ-1) with a charge-coupled device (CCD) sensor. CHRIS/PROBA = Compact High-Resolution Imaging Spectrometer/Project for On-Board Autonomy; MODIS = moderate-resolution imaging spectroradiometer; AWiFS = Advanced Wide-Field Sensor

Table 1. Summary of LAI retrieval based on LUTs generated from the PROSAIL model.

LUTs	Platforms	Sensors	References
Reflectance-LUTs	UAV/Airborne	Multispectral camera	[12,15]
		Hyperspectral camera	[11,13]
	Satellite	ZY-3 MUX, GF-1 WFV, HJ-1 CCD	[18]
		Landsat	[19,20,21]
		Sentinel-2	[22,23,24]
		RapidEye	[22]
		CHRIS/PROBA	[25]
		MODIS	[20,26]
		AWiFS	[27]
	Ground-based	Hyperspectral spectrometer	[28,29,30]
VI-LUTs	satellite	CHRIS/PROBA	[31]
		RapidEye	[17]
		Landsat	[32]

LAI = Leaf area index; UAV = unmanned aerial vehicle; Reflectance-LUTs = reflectance-based lookup tables; VI-LUTs = Vegetation index-based lookup tables; ZY-3 MUX = the resources satellite No. 3 (ZY-3) with a multi-spectral camera (MUX); GF-1 WFV = the satellite GaoFen No. 1 (GF-1) with a wide field of view (WFV) camera; HJ-1 CCD = the environment satellite (HJ-1) with a charge-coupled device (CCD) sensor. CHRIS/PROBA = Compact High-Resolution Imaging Spectrometer/Project for On-Board Autonomy; MODIS = moderate-resolution imaging spectroradiometer; AWiFS = Advanced Wide-Field Sensor

The radiometric data acquired from various sensors have different spectral and spatial characteristics because of their sensor-dependence [3]. Previous studies have already demonstrated that normalized difference vegetation index (NDVI) is significantly affected by differences in bandwidth and spatial resolutions [33,34]. Owing to weight reduction and downsizing, spectral cameras onboard UAVs have technical differences in the central wavelengths and bandwidths compared with cameras onboard ground and satellite platforms. Most VIs were developed based on the satellite (e.g., enhanced vegetation index without a blue band (EVI2) [35], atmospherically resistant vegetation index (ARVI) [36]) and ground-based remote-sensing data (e.g., soil adjusted vegetation index (SAVI), modified soil adjusted vegetation index (MSAVI) [37]). Thus the calculation of these VIs using field-scale UAV images might cause uncertainties for crop parameter estimation. Red (R) and near-infrared (NIR) bands have been widely used for calculations of VIs because the spectral reflectance at R and NIR bands is sensitive to vegetation phenotypic traits (e.g., LAI and chlorophyll content) [38,39,40,41]. Previous studies have attempted to optimize the construction of VIs by evaluating the central wavelengths and bandwidths of R and NIR bands [39]. Therefore, this study followed this method to improve the performance of UAV-based VIs on LAI estimation at the field scale.

The objectives of this study were to (1) search for a feasible LUT strategy for field-scale LAI retrieval using multispectral UAV images, and (2) improve the easiness to use, robustness, and accuracy of this method so that it can provide precise agronomic information at field scale. To that purpose, the performance of reflectance-LUT and VI-LUT strategies generated by the PROSAIL model for LAI retrieval were investigated using multispectral UAV images; the central wavelengths and bandwidths of the red and NIR bands for VIs calculations on LAI retrieval were evaluated and further optimized using UAV hyperspectral images.

2. Methodology

2.1. Study Area and Long-Term Experimental Plots

UAV observations were conducted at the Yucheng Comprehensive Experiment Station (YCES) of the Chinese Academy of Sciences (36.83° N, 116.57° E), which is located along the north side of the Lower Yellow River (Figure 1). The study area has a warm temperate and a semi-humid monsoon climate with an approximate annual mean temperature of 13.4 °C and an average annual precipitation of 576.70 mm, concentrated between July and September. The YCES is dominated by a typical cropping system of winter wheat (October–June) and summer maize (July–September).

We selected six types of long-term experimental plots with different tillage treatments, soil water content, and nutrient levels for experiments (Figure 1 and

Table 2. Treatments of the six experimental fields.

Field (Plot Number, Plot Size)	Treatments 1 (Plot Number)
Field A (16, 2×2 m2)	CT-RR 2 (3), NT-RR (3), CT-S (3), NT-S (4), CT-1/2S (3)
Field B (25, 5×6 m2)	CK (4), NP (4), NK (4), PK (4), NPK (7), NPK-S (2)
Field C (6, 10×10 m2)	CT (6)
Field D (32, 5×10 m2)	N0-60%fc (3), N70-60%fc (3), N140-60%fc (3), N210-60%fc (3) N280-60%fc (4), N0-80%fc (3), N70-80%fc (3), N140-80%fc (3) N210-80%fc(3), N280-80%fc(4)
Field E (32, 5×10 m2)	0%fc (8), 40%fc (8), 60%fc (8), 80%fc (8)
Field F (18, 8×40 m2)	NT-S-F1 (3), NT-S-F2 (3), NT-RR-F1 (3), NT-RR-F2 (3), CT-S-F1 (3), CT-S-F2 (3)

¹ Tillage treatments: see Table 3; ² The symbol “-” represents different treatments are implemented together in the plot.

). All the experimental plots had been operational for >10 years. Owing to the obvious differences in crop traits caused by the existing gradients of agricultural treatments, the crop LAI dynamic range was large. Thus, the study area was perfect for UAV remote-sensing experiments.

2.2. Ground Measurements and UAV Flight Missions

2.2.1. Ground Measurements of LAI

Leaves of wheat were collected via destructive sampling for LAI measurements in the laboratory (Figure 2). The row spacing of all experimental plots was 20 cm. In order to ensure the reasonability and typicality of the study, wheat with relatively uniform canopies outside the border areas of each plot was selected for the experiment. A row of wheat with a length of 0.5 m was cut for LAI measurements. The leaf area was measured using an LI-3000C Leaf Area Meter (Li-COR Biosciences, Lincoln, NE, USA). The average LAI of each plot was calculated according to

LAI (m²/m²) = TLA/(a × b)

(1)

where TLA (m²) represents the total leaf area of crops within each plot, and a (m) and b (m) represent the width and length of each plot, respectively. Because wheat is row crops, the TLA was calculated according to

TLA = 2 × LA × nrows × b

(2)

where LA (m²) represents the measured leaf area of crops harvested within each plot, nrows represents the number of rows for each plot. With UAV flights, three LAI datasets were collected for evaluation of the LAI retrieval accuracy (

Table 4. Statistical summary of the in situ measured LAI (m²/m²).

Cameras	Date	Field (Samples)	Min	Mean	Max	Standard Deviation	CV (%) *
Multi-	15 May 2018	A (4)	1.65	1.78	1.85	0.10	5.40
		B (25)	0.51	1.27	2.30	0.59	46.34
		C (6)	2.08	2.59	3.05	0.37	14.18
		D (32)	0.31	1.83	3.15	0.77	42.20
		E (32)	0.63	1.37	1.85	0.30	22.16
		F (8)	1.78	2.06	2.50	0.25	11.92
		Total (107)	0.31	1.62	3.15	0.65	40.03
	16 May 2019	B (25)	0.21	1.25	2.56	0.79	63.30
		D (32)	0.46	1.51	2.93	0.78	51.73
		Total (57)	0.21	1.40	2.93	0.81	57.70
Hyper-	15 May 2018	B (25)	0.51	1.27	2.30	0.59	46.34
		D (32)	0.31	1.83	3.15	0.77	42.20
		E (32)	0.63	1.37	1.85	0.30	22.16
		Total (89)	0.31	1.51	3.15	0.63	41.85

* CV represents the coefficient of variation.

). Each dataset exhibited a high coefficient of variation (CV) owing to the wide range of agricultural treatments, which contributed to a good model fit between the measured and simulated LAI. Because the flight coverage of the multispectral UAV observations (all six fields) differed from that of hyperspectral UAV observations (see Figure 1 for the two red rectangles), the number of ground samples for the LAI estimation was different for 15 May 2018.

2.2.2. UAV Flight Missions and Data Pre-Processing

On 15 May 2018, and 16 May 2019, two multispectral UAV flight observations were conducted. On 15 May 2018, one hyperspectral UAV flight observation was conducted. For the UAV missions, two types of UAV platforms (see Figure 3 for the DJ M100 four-rotator UAV (SZ DJI Technology Co., Shenzhen, Guangzhou, China) and DJ M600 Pro six-rotator UAV (SZ DJI Technology Co., Shenzhen, Guangzhou, China) equipped with two sensors (

Table 5. Specifications of the two cameras used in this study.

Sensor	Spectral Channels (Central Wavelength/Spectral Ranges)	Spatial Resolution	Spectral Resolution (nm)
Cubert S185	450–950 nm	1 cm	4 nm
RedEdge-M	B475, G560, R668, E717, NIR840	4 cm	B (20), G (20), R (10), E (10), NIR (40)

)—Micasense RedEdge-M multispectral camera (MicaSense, Seattle, WA, USA) and Cubert S185 hyperspectral camera (Cubert GmbH, Ulm, Germany)—were used. The two multispectral flights covered all six fields; the hyperspectral flight covered Field B, D, and E (see Figure 1 for the two red rectangles).

The UAV flight time of the three flight missions was from 10:00 to 14:00. For radiation correction, the reflectance of a spectral panel was collected during each flight. The radiation correction, image mosaicking, and orthography of the UAV images were conducted using Agisoft Photoscan (Agisoft LLC, St. Petersburg, Russia) and Pix4D Mapper 3.1.22 (Pix4D, S.A., Lausanne, Switzerland). ENVI 5.1 (ESRI, RedLands, CA, USA) and Python 2.7 were used for further data analyses and programming, respectively.

2.3. Retrieving LAI from UAV Data Using PROSAIL Model

The PROSAIL model, along with the UAV remote-sensing data, was used to retrieve wheat LAI [42]; it is derived from the combination of the PROSPECT blade model and the SAIL canopy structure model [43]. The PROSPECT model simulates the optical properties of leaves, from 400 to 2500 nm, with four inputs: leaf structure parameter (N), chlorophyll content (Chl), leaf water mass per area (LMA), and blade equivalent thickness (EWT). The SAILH model is a radiative transfer model on the canopy scale, in which vegetation is treated as a mixed medium with an assumption that the blade azimuth distribution is uniform [44]. In total, fourteen parameters (N, Chl, EWT, LMA, LAI, leaf carotenoid content (caro), brown pigment content, soil brightness parameter (psoil), hot-spot size parameter (hot spot), solar zenith and azimuth angles, view zenith and azimuth angles, and average leaf angle (ALA)) are required for running the PROSAIL model, most of which are difficult to obtain (

Table 6. Ranges of input variables for the PROSAIL model.

Variable	Abbr.	Unit	Value (LUT) *	Range (EFAST) **
Leaf structure parameter	N	Unitless	1.5	1–2
Leaf chlorophyll content	Chl	μg·cm−2	20–70 (step = 0.2)	20–70
Leaf carotenoid content	caro	μg·cm−2	10	3–30
Brown pigment content	-	arbitrary units	0	0
Blade equivalent thickness	EWT	cm	0.01	0.005–0.03
Leaf water mass per area	LMA	g·cm−2	0.005	0.004–0.007
Soil brightness parameter	psoil	Unitless	0.1	0.01–0.3
Leaf area index	LAI	m2 m−2	0.1–6 (step = 0.01)	0.1–6
Hot-spot size parameter	hot spot	m m−1	0.2	0.05-1
Solar zenith angle	-	degrees	20	20
Solar azimuth angle	-	degrees	185	185
View zenith angle	-	degrees	0	0
View azimuth angle	-	degrees	0	0
Average leaf angle	ALA	degrees	70	30–70

* Value (LUT) is used for generating lookup tables in Section 2.3.2; ** Range (Fourier amplitude sensitivity test (EFAST)) was used for running the global sensitivity analyses of PROSAIL model in Section 2.3.1.

) (e.g., ALA, caro, and hot spot).

High-dimensional LUTs are often generated to retrieve LAI because most PROSAIL inputs are difficult to reach in practice [11,24]. The computation processes for high-dimensional LUTs are slow and complex, which may hinder the application of LUTs to parameter inversion. Sensitivity analysis of models is an alternative way to solve this issue. It can identify the sensitive model inputs and evaluate their sensitivity levels; the changes of these sensitive inputs in a certain range will lead to obvious variations of model outputs. Thus some unessential PROSAIL variables can be set as constants so that the dimension of the LUTs decreases [12]. Sensitivity analysis includes local and global analyses. The global analysis can reveal the effects of each input and interaction among inputs on model outputs. In this study, the extended Fourier amplitude sensitivity test (EFAST) was used for global sensitivity analysis. The EFAST is a quantitative method based on variance [45]. The variance produced by the change of model outputs reveals the sensitive values or levels of model inputs [12].

LAI retrieval from the UAV data comprised three steps (Figure 4). Firstly, select appropriate VIs for LAI retrieval based on the (EFAST) global sensitivity analyses of the PROSAIL model; secondly, generate reflectance-LUTs and VI-LUTs by running the PROSAIL model; thirdly, retrieve LAI through cost functions (see Section 2.3.3 for the formula of cost functions) based on the LUTs. In step 3, both of the multispectral and hyperspectral UAV datasets were used. The multispectral datasets were used to investigate the robustness of LUT strategies for LAI retrieval (step 3.1); the hyperspectral dataset was used to improve the LAI retrieval robustness and accuracy (step 3.2).

2.3.1. Selecting Optimal VIs for LAI Retrieval (Global Sensitivity Analysis)

The global sensitivity analyses were conducted to select VIs for LAI retrieval. First, 5000 model input datasets (Range (EFAST), Table 6) with uniform distributions were generated using the Simlab 2.2.1 software (JRC, Italy). These generated datasets were the simulated inputs of the PROSAIL model. Second, the PROSAIL model was run in a forward mode to obtain model outputs (spectral reflectance) using the PyProSAIL package of Python 2.7 (http://teledetection.ipgp.jussieu.fr/prosail/, accessed on 24 June 2019). Then VIs (Table 7) were calculated using the PROSAIL outputs—spectral reflectance. Thirdly, datasets of the PROSAIL inputs and outputs were used to conduct the EFAST global sensitivity analysis using Simlab 2.2.1.

The reflectance values of the blue (B), green (G), red (R), red-edge (E), and NIR bands were denoted as ρ_B, ρ_G, ρ_R, ρ_E, and ρ_NIR, respectively. We selected some VIs that were used for the estimation of crop traits in previous studies (Table 7, where ρ₁ = ρ_NIR). These VIs were divided into three groups: R-VIs (calculated using ρ_R and ρ_NIR), G-VIs (calculated using ρ_G and ρ_NIR), and E-VIs (calculated using ρ_E and ρ_NIR). In this study, variables of VIs were changed; thus VIs with the same formula but different variables were denoted as m-VIs (e.g., modified atmospherically resistant vegetation index (m-ARVI) represents R-ARVI (ρ₂ = ρ_R), G-ARVI (ρ₂ = ρ_G), and E-ARVI (ρ₂ = ρ_E) in the R-VI, G-VI, and E-VI groups, respectively).

Table 7. Bands and VIs used in this study.

m-VIs	Equations	R-VIs	G-VIs	E-VIs	Ref.
		ρ2 = ρR	ρ2 = ρE	ρ2 = ρE
m-ARVI	(ρ1 − 2ρ2 + ρB)/(ρ1 + 2ρ2 − ρB)	√	√	√	[3]
m-EVI	2.5 ((ρ1 − ρ2)/(ρ1 + 6ρ2 − 7.5ρB + 1))	√	√	√	[46]
m-WDRVI	(0.12ρ1 − ρ2)/(0.12ρ1 + ρ2)	√	√	√	[46]
m-MSR	$({\rho }_1/{\rho }_2-1)/\sqrt{\left({\rho }_1/{\rho }_2+1\right)}$	√	√	√	[47,48]
m-MSAVI2	$0.5[\left(2{\rho }_1+1\right)-\sqrt{{\left(2{\rho }_1+1\right)}^2-8\left({\rho }_1-{\rho }_2\right)}]$	√	√	√	[38]
m-MCARI2	$1.5[2.5\left({\rho }_1-{\rho }_2\right)-1.3\left({\rho }_1-{\rho }_G\right)]/\sqrt{[{\left(2{\rho }_1+1\right)}^2-\left(6{\rho }_1-5{\rho }_2\right)-0.5]}$	√	√	√	[48]
m-MTVI1	1.2(1.2(ρ1 − ρG) − 2.5(ρ2 − ρG))	√		√	[48]
m-TVI	0.5(120(ρ1 − ρG)) − 200(ρ2 − ρG)	√		√	[48]
m-NDVI	(ρ1 − ρ2)/(ρ1 + ρ2)	√	√	√	[28]
m-OSAVI	1.16(ρ1 − ρ2)/(ρ1 + ρ2 + 0.16)	√	√	√	[46]
m-RVI	ρ1/ρ2	√	√	√	[3]
m-SAVI	(1 + 0.5)((ρ1 − ρ2)/(ρ1 + ρ2 + 0.5))	√	√	√	[3]
m-NRI	ρ2/(ρ2 + ρE + ρ1)	√	√		[49]
m-EVI2	2.5(ρ1 − ρ2)/(ρ1 + 2.4ρ2 + 1)	√	√	√	[35,46]

R = Red, G = Green, E = red-edge; m-ARVI = modified atmospherically resistant vegetation index; m-EVI = modified enhanced vegetation index; m-WDRVI = modified wide dynamic range vegetation index; m-MSR = modified simple ratio; m-MSAVI2 = modified secondary optimized soil adjusted vegetation index; m-MCARI2 = modified improved chlorophyll absorption ratio index; m-MTVI1 = modified improved triangular vegetation index; m-TVI = modified triangular vegetation index; m-NDVI = modified normalized difference vegetation index; m-OSAVI = modified optimized soil adjusted vegetation index; m-RVI = modified ratio index; m-SAVI = modified soil adjusted vegetation index; m-NRI = modified nitrogen ratio index; m-EVI2 = modified enhanced vegetation index without a blue band.

Table 7. Bands and VIs used in this study.

m-VIs	Equations	R-VIs	G-VIs	E-VIs	Ref.
		ρ2 = ρR	ρ2 = ρE	ρ2 = ρE
m-ARVI	(ρ1 − 2ρ2 + ρB)/(ρ1 + 2ρ2 − ρB)	√	√	√	[3]
m-EVI	2.5 ((ρ1 − ρ2)/(ρ1 + 6ρ2 − 7.5ρB + 1))	√	√	√	[46]
m-WDRVI	(0.12ρ1 − ρ2)/(0.12ρ1 + ρ2)	√	√	√	[46]
m-MSR	$({\rho }_1/{\rho }_2-1)/\sqrt{\left({\rho }_1/{\rho }_2+1\right)}$	√	√	√	[47,48]
m-MSAVI2	$0.5[\left(2{\rho }_1+1\right)-\sqrt{{\left(2{\rho }_1+1\right)}^2-8\left({\rho }_1-{\rho }_2\right)}]$	√	√	√	[38]
m-MCARI2	$1.5[2.5\left({\rho }_1-{\rho }_2\right)-1.3\left({\rho }_1-{\rho }_G\right)]/\sqrt{[{\left(2{\rho }_1+1\right)}^2-\left(6{\rho }_1-5{\rho }_2\right)-0.5]}$	√	√	√	[48]
m-MTVI1	1.2(1.2(ρ1 − ρG) − 2.5(ρ2 − ρG))	√		√	[48]
m-TVI	0.5(120(ρ1 − ρG)) − 200(ρ2 − ρG)	√		√	[48]
m-NDVI	(ρ1 − ρ2)/(ρ1 + ρ2)	√	√	√	[28]
m-OSAVI	1.16(ρ1 − ρ2)/(ρ1 + ρ2 + 0.16)	√	√	√	[46]
m-RVI	ρ1/ρ2	√	√	√	[3]
m-SAVI	(1 + 0.5)((ρ1 − ρ2)/(ρ1 + ρ2 + 0.5))	√	√	√	[3]
m-NRI	ρ2/(ρ2 + ρE + ρ1)	√	√		[49]
m-EVI2	2.5(ρ1 − ρ2)/(ρ1 + 2.4ρ2 + 1)	√	√	√	[35,46]

R = Red, G = Green, E = red-edge; m-ARVI = modified atmospherically resistant vegetation index; m-EVI = modified enhanced vegetation index; m-WDRVI = modified wide dynamic range vegetation index; m-MSR = modified simple ratio; m-MSAVI2 = modified secondary optimized soil adjusted vegetation index; m-MCARI2 = modified improved chlorophyll absorption ratio index; m-MTVI1 = modified improved triangular vegetation index; m-TVI = modified triangular vegetation index; m-NDVI = modified normalized difference vegetation index; m-OSAVI = modified optimized soil adjusted vegetation index; m-RVI = modified ratio index; m-SAVI = modified soil adjusted vegetation index; m-NRI = modified nitrogen ratio index; m-EVI2 = modified enhanced vegetation index without a blue band.

2.3.2. Generating Reflectance-LUTs and VI-LUTs

Since sensitivity analysis can identify sensitive inputs of a model and evaluate their sensitivity levels, some unessential inputs can be set as constant values for decreasing the dimension of LUTs. For simplifying the LAI retrieval to fulfill the requirement of practical application, we generated two-dimensional (2D) LUTs for LAI retrieval, as previous studies have done [12], which set chlorophyll content and LAI as variables while other parameters in PROSAIL model as constants.

The ranges of the model inputs were determined for wheat according to ground measurements, previous studies [23], and the LOPEX’93 database (http://opticleaf.ipgp.fr/index.php?page=database, accessed on 29 November 2018) (Value (LUT), Table 6). We generated two types of LUTs: the reflectance-LUTs and VI-LUTs, which retrieved the LAI through spectral reflectance and VIs, respectively.

2.3.3. Retrieving LAI through Cost Functions

The cost function was used to find the optimal LAI estimate

$$f_{\mathrm{RMSE}}=\sqrt{\frac{1}{p}{{\displaystyle \sum }}_{i=1}^p{\left(R_m-R_s\right)}^2}$$

(3)

where p represents the number of input variables (the number of bands for reflectance-LUTs/p = 1 for VI-LUTs). R_m represents the measured spectral reflectance/VIs, and R_s represents the simulated reflectance/VIs from the PROSAIL model. When the cost-function value reaches the minimum in a given reflectance-LUT/VI-LUT, the simulated LAI is considered to be the optimal estimate.

In this section, both of the multispectral and hyperspectral UAV data were used for LAI retrieval; the multispectral data were used for investigating the robustness of LAI retrieval strategies. However, the hyperspectral data were further used for improving LAI retrieval accuracy and robustness. For clear statements, the optimization of LAI retrieval is explained in detail in Section 2.3.4.

2.3.4. Optimizing LAI Retrieval Using Hyperspectral Datasets

The hyperspectral UAV data were used for optimizing LAI retrieval. The central wavelengths and bandwidths of the R and NIR bands for LAI retrieval were assessed to optimize the construction of four UAV-based VIs (m-NDVI, m-MCARI2, m-ARVI, and m-NRI selected by global sensitivity analysis; see Section 3.1); the central wavelengths of the B, G, and E bands were set to constant values of 475, 560, and 717 nm, respectively, according to the specifications of the Micasense RedEdge-M (Table 5).

For evaluating the effects of the central wavelengths and bandwidth on LAI retrieval, four simulated datasets were generated as

Table 8. Four UAV datasets used for LAI retrieval.

Dataset	Value of Central Wavelengths (Bandwidth) (nm)
	B	G	R	E	NIR
Dataset 1	475/20	560/20	668/10	717/10	840/40
Dataset 2	475/20	560/20	672 or 612/10 *	717/10	752 or 756/10 **
Dataset 3	475/4	560/4	668/4	717/4	840/4
Dataset 4	475/4	560/4	672 or 612/10 *	717/4	752 or 756/10 **

* The central wavelengths of R bands were 672 nm for m-NDVI, m-NRI, and m-ARVI, and 612 nm for m-MCARI2. Please see Table 10 for the selection of optimal central wavelengths. ** The central wavelengths of NIR bands were 752 nm for m-NDVI, m-NRI, and m-ARVI, and 756 nm for m-MCARI2. Please see Table 10 for the selection of optimal central wavelengths.

shows. The bandwidths of Datasets 1 and 2 were broad; those of Datasets 3 and 4 were narrow (4 nm). The central wavelengths of R and NIR bands were optimized in Datasets 2 and 4.

The optimal central wavelengths of the NIR and R bands were determined via three steps. First, an autocorrelation analysis was performed between any two hyperspectral bands among 88 bands (600–950 nm, Cubert S185). We set p > 0.05 to find uncorrelated bands since the p value of correlated variables was < 0.05. Second, four selected m-VIs (m-NDVI, m-MCARI2, m-ARVI, and m-NRI; see Section 3.1) were calculated from the combinations of reflectance data determined above. Third, the optimal central wavelengths were determined by the highest Pearson correlation coefficient (r values) between the calculated VIs and the measured LAI.

Then, for evaluating the effects of the bandwidth on the LAI retrieval, the narrow bandwidths of the hyperspectral data were resampled into broad bandwidths (Table 8). Because the spectral resolution of the hyperspectral data (Cubert S185) was relatively coarse (4 nm), we only compared the LAI retrieval accuracy of VIs with a narrow bandwidth (4 nm) and broad bandwidths (equal to those of the RedEdge-M camera; see Table 5).

2.4. Statistical Analysis

Python 2.7 was used for statistical analyses, including four statistical indicators for evaluation of the LAI retrieval accuracy: the Pearson correlation coefficient (r), coefficient of determination (R²), root-mean-square error (RMSE), and mean relative error (MRE). The equations for calculating the R², RMSE, and MRE were calculated according to

$$R^2={\{[{{\displaystyle \sum }}_{i=1}^n\left(M_i-\overline{M}\right)\left(E_i-\overline{E}\right)]/[\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(M_i-\overline{M}\right)}^2}\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(E_i-\overline{E}\right)}^2}]\}}^2$$

(4)

$$\mathrm{RMSE}=\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(M_i-E_i\right)}^2/n}$$

(5)

$$\mathrm{MRE}=[{\displaystyle \sum }_{i=1}^n\left(\left|M_i-E_i\right|/M_i\right)]/n$$

(6)

where i represents the sequence number of the array, M_i represents the LAI value of the measured LAI array, E_i represents the LAI value of the estimated LAI array, $\overline{M}$ represents the average value of the measured LAI array, and $\overline{E}$ represents the average value of the estimated LAI array.

3. Results

3.1. Optimal VIs Selected through Global Sensitivity Analyses

Results of the global sensitivity analyses indicate that the B, R, and NIR bands were more sensitive to LAI, whereas the G and E bands were more sensitive to Chl (Figure 5a). Additionally, E and NIR bands were quite sensitive to the average leaf angle (ALA), with values of 0.26 and 0.19, respectively; B and R were quite sensitive to soil-brightness parameter (psoil), with values of 0.35 and 0.27, respectively; G and E were partially sensitive to leaf-structure parameter (N), with values of 0.17 and 0.11, respectively. Compared with the multispectral bands, the average sensitivity values of the VIs to psoil and N decreased to 0.02 and 0.03, respectively. Thus, the LAI, Chl, and ALA were the three most sensitive PROSAIL inputs to the 39 VIs, with average sensitivity values of 0.72, 0.21, and 0.10, respectively. The average sensitivity values of the R-VIs to the LAI, Chl, and ALA were 0.86, 0.03, and 0.13, respectively. Those of the E-VIs were 0.60, 0.39, and 0.06, respectively, and those of the G-VIs were 0.71, 0.22, and 0.10, respectively. The foregoing results indicate that the R-VIs were the most sensitive to the LAI, the least sensitive to the Chl. Therefore, R-VIs were selected for LAI retrieval.

Within the R-VI group, the total sensitivity values of the following six R-VIs to the LAI reached 0.90: nitrogen ratio index in the R-VI group (R-NRI, 0.93), normalized difference vegetation index in the R-VI group (R-NDVI, 0.93), wide dynamic range vegetation index in the R-VI group (R-WDRVI, 0.91), modified improved chlorophyll absorption ratio index in the R-VI group (R-MCARI2, 0.91), atmospherically resistant vegetation index in the R-VI group (R-ARVI, 0.91), and optimized soil adjusted vegetation index in the R-VI group (R-OSAVI, 0.90). R-NDVI, R-WDRVI, and R-OSAVI were calculated using only the R and NIR bands, while R-NDVI was more commonly used and more sensitive to the LAI than R-WDRVI and R-OSAVI. Therefore, R-NDVI, R-NRI, R-MCARI2, and R-ARVI were selected for LAI retrieval. For simplicity, R-NDVI, R-NRI, R-MCARI2, and R-ARVI are denoted as NDVI, NRI, MCARI2, and ARVI, respectively.

3.2. LAI Retrieval Based on Two LUT Strategies Using Multispectral UAV Data

For assessment of the LAI retrieval performance of two LUT strategies, two-year UAV multispectral observations and were used. For reflectance-LUT, the R² values were 0.42 and close to zero in 2018 and 2019, respectively, remarkably lower than those for the VIs-LUT (R² > 0.74) (

Table 9. R², RMSE, and MRE for linear regression between measured LAI and simulated LAI based on reflectance-LUTs/VI-LUTs using multispectral UAV data.

Year	Reflectance-LUTs				VI-LUTs
	Bands	R2	RMSE	MRE	VIs	R2	RMSE	MRE
2018 n = 107	R,NIR	0.42	0.94	0.70	NDVI	0.76	0.44	0.25
	R,NIR,B	0.42	0.94	0.70	ARVI	0.74	0.51	0.30
	R,NIR,G	0.42	0.94	0.70	MCARI2	0.75	0.38	0.22
	R,NIR, E	0.42	0.94	0.70	NRI	0.75	0.46	0.25
2019 n = 57	R,NIR	0.01	2.23	2.71	NDVI	0.78	0.38	0.27
	R,NIR,B	0.01	2.23	2.71	ARVI	0.74	0.47	0.23
	R,NIR,G	0.01	2.23	2.71	MCARI2	0.83	0.33	0.30
	R,NIR, E	0.01	2.23	2.71	NRI	0.74	0.43	0.31

R² represents the coefficient of determination for a linear regression model between the measured LAI and simulated LAI. RMSE = root-mean-square error, MRE = mean relative error.

). The RMSE and MRE for the reflectance-LUT (RMSE ≥ 0.94, MRE ≥ 0.70 in 2018; RMSE ≥ 2.23, MRE ≥ 2.71 in 2019) were higher than those for the VIs-LUT (RMSE ≤ 0.51, MRE ≤ 0.30 in 2018; RMSE ≤ 0.47, MRE ≤ 0.31 in 2019). The foregoing results indicate that the VI-LUT strategy was more robust and accurate than the reflectance-LUT strategy, which offers an alternative method for LAI retrieval.

The LAI retrieval from four VI-LUTs corresponded to ground-based measurements (Figure 6 and Table 9). The R², RMSE, and MRE values for LAI retrieval ranged between 0.74–0.83, 0.33–0.51, and 0.22–0.31, respectively. Overall, the LAI estimation accuracy was similar among the four VI-LUTs, although the NDVI-LUT and MCARI2-LUT performed slightly better. In 2018, the R² values of the MCARI2-LUT and NDVI-LUT were close to 0.75, and the RMSE and MRE of MCARI2-LUT were 0.38 and 0.22, respectively (the lowest among the four VI-LUTs), which was consistent with the comparison between the regression line and the 1:1 line (Figure 6c). In 2019, the R² value of the MCARI2-LUT was 0.83, higher than the other three VI-LUTs. The RMSE and MRE values of MCARI2-LUT were 0.33 and 0.30, respectively, which were the lowest among the four VI-LUTs. For ARVI-LUT and NRI-LUT, LAI values were slightly underestimated both in 2018 and in 2019. The foregoing results indicate that the VI-LUTs had different characteristics of LAI estimation.

3.3. Optimization of VI-LUTs for LAI Retrieval Using Hyperspectral UAV Data

3.3.1. Optimization of Central Wavelengths for VI Calculation

Figure 7a shows the results of the autocorrelation analyses (Pearson correlation analysis between any two bands). Figure 7b is the uncorrelated band ranges under the conditions of r < 0.20 and p > 0.05 (n = 89, see Table 4). As the wavelength that located in NIR ranges increases, the signal-noise ratio gradually becomes low and the reflectance will be more easily affected by the absorption of water [50]. Thus, the band ranges over 850 nm were excluded in this study. We selected two appropriate regions: 740–760 nm (ρ₁ bands) and 600–723 nm (ρ₂ bands), along with 740–850 nm (ρ₁ bands) and 703–725 nm (ρ₂ bands).

Then, m-NDVI, m-MCARI2, m-ARVI, and m-NRI were calculated using the combinations of ρ₁ bands and ρ₂ bands within the above-determined ranges. The central wavelengths of the ρ₁ and ρ₂ bands were refined via Pearson correlation analysis between the calculated VIs and the measured LAI (Figure 8). The specific bands of ρ₁ and ρ₂ were determined for four VIs by identifying the highest r value from the correlation plots in Figure 8. As shown in

Table 10. Specific optimal bands of ρ₁ and ρ₂ for the four m-VI calculations.

Crops	VIs	ρ1 Bands (nm)	ρ2 Bands (nm)	r Value
Wheat	m-NDVI	752	672	0.86 **
	m-NRI	752	672	−0.86 **
	m-MCARI2	756	612	0.87 **
	m-ARVI	752	672	0.86 **

** indicates p < 0.01.

, m-NDVI, m-MCARI2, and m-ARVI exhibited significantly positive correlations with the measured LAI (r > 0.86, p < 0.01, n = 89), while m-NRI was significantly negatively correlated with the measured LAI (r_m-NRI = −0.86, p < 0.01, n = 89). Correspondingly, the optimal central wavelengths of m-NDVI, m-ARVI, and m-NRI were all located at 752 nm (ρ₁ band) and 672 nm (ρ₂ band), while the values for m-MCARI2 were located at 756 nm (ρ₁ band) and 612 nm (ρ₂ band).

3.3.2. LAI Retrieval Based on Two LUT Strategies Using Hyperspectral UAV Data

A comparison between the measured and simulated LAI indicates that the VI-LUT strategy with UAV hyperspectral observations was more robust than the reflectance-LUT strategy (

Table 11. R², RMSE, and MRE for linear regression between simulated LAI and estimated LAI based on reflectance-LUTs/VI-LUTs using hyperspectral UAV data in 2018 (n = 89).

Reflectance-LUTs				VIs-LUTs
Bands	R2	RMSE	MRE	VIs	R2	RMSE	MRE
R, NIR	0.27	2.58	2.22	NDVI	0.80	0.55	0.31
R, NIR, B	0.27	2.58	2.22	ARVI	0.76	0.66	0.40
R, NIR, G	0.27	2.58	2.22	MCARI2	0.82	0.37	0.26
R, NIR, E	0.32	2.60	2.24	NRI	0.75	0.58	0.32

). The R² values for the reflectance-LUTs (≤0.32) were significantly lower than those for the VIs-LUTs (≥0.75), and the RMSE and MRE values for the reflectance-LUTs (RMSE > 2.57, MRE ≥ 2.22) were higher than those for the VIs-LUTs (RMSE ≤ 0.66, MRE ≤ 0.40). Furthermore, among the four proposed VIs, the MCARI2-LUT was the most robust and accurate strategy for wheat LAI retrieval (R² = 0.82, RMSE = 0.37, and MRE = 0.26).

3.3.3. Evaluation of Optimized VI-LUTs Using Hyperspectral Data for LAI Retrieval

Firstly, we compared datasets with the same bandwidths to reveal the effects of central wavelengths on LAI retrieval. For datasets with a broad bandwidth (Datasets 1 and 2, the dots; see Figure 9), the LAI retrieval robustness for the MCARI2-LUT was improved slightly after optimizing the central wavelengths of the VIs, as its RMSE decreased by 0.22. However, those of ARVI-LUT did not. For datasets with a narrow bandwidth (Datasets 3 and 4, the crosses), after optimizing central wavelengths, R² increased obviously by 0.12, 0.10, 0.12, and 0.09 for the NDVI-LUT, ARVI-LUT, MCARI2-LUT, and NRI-LUT, respectively. RMSE decreased by 0.04 and 0.06, and MRE decreased by 0.03 and 0.06 for NDVI-LUT and MCARI2-LUT, respectively. Although RMSE of ARVI-LUT and NRI-LUT slightly increased by 0.02 and 0.01, respectively, overall, optimizing the central wavelengths of datasets with narrow bandwidths could improve the robustness of VI-LUTs for LAI retrieval.

Secondly, we compared datasets with the same central wavelengths to show the effects of bandwidth on LAI retrieval. For datasets with the optimized central wavelengths (Datasets 2 and 4, in red, see Figure 9), overall, the LAI retrieval for datasets with narrow bandwidths was more robust than datasets with broad bandwidths. However, for VIs without optimization of the central wavelengths (Datasets 1 and 3, in black; see Figure 9), overall, the effects of bandwidth on LAI retrieval remained unclear according to the results of this study. These foregoing results indicate that both of the central wavelength and bandwidth affected LAI retrieval.

4. Discussion

Our research indicates that low-dimensional VI-LUTs along with UAV multispectral images were feasible for wheat LAI retrieval, and more robust than reflectance-LUTs. The MCARI2-LUT and NDVI-LUT are recommended for LAI retrieval owing to their better performance. Furthermore, the central wavelengths and bandwidths of a commonly used multispectral camera (Micasense RedEdge-M), can be modified for specific applications. Regarding the LAI retrieval accuracy with two-year multispectral UAV images, the R² values of the VI-LUT strategy were >0.74. The R² values of the reflectance-LUT strategy were 0.42 in 2018 and close to zero in 2019. In wheat LAI retrieval with four simulated datasets, a dataset with central wavelengths of 756 nm for the NIR band and 612 nm for the R band and a narrow bandwidth (~4 nm) performed the most robustly; the R² value of the MCARI2-LUT increased by 0.06, and the RMSE and MRE decreased by 0.18 and 0.10, respectively.

4.1. Analyses of LAI Retrieval Performance for Reflectance-LUTs and VI-LUTs

Results of the global sensitivity analyses for PROSAIL model indicate that the LAI was the dominant factor affecting B, R, and NIR bands and VIs, in agreement with previous studies [43,51]. The Chl and ALA were two additional variables that the spectral bands and VIs were sensitive to (see Figure 5). In this study, the VIs were divided into three groups to identify the optimal VI group for LAI retrieval. The results indicate that the R-VIs were the most sensitive to the LAI and could reduce the influences of other PROSAIL variables. The E-VIs and G-VIs were sensitive to the Chl but less influenced by the ALA. This study selected the R-VIs for LAI retrieval. However, E-VIs have also been widely applied in LAI retrieval owing to their capabilities of reducing the impact of the ALA of different crops [17]. But the interactions between crop biophysical (e.g., LAI, ALA) and biochemical parameters (e.g., Chl) may introduce uncertainties in parameter retrieval [28]. Additionally, the red-edge bandwidth is narrow, with a low signal-to-noise ratio and uncertainties in its spectral response [51]. Therefore, in UAV remote-sensing, replacing the traditional red bands (ρ₂ = ρ_R) with red-edge bands (ρ₂ = ρ_E) for VIs calculation on LAI retrieval must be considered in future studies.

The results based on two-year multispectral data indicate that the VI-LUT strategy was able to accurately estimate the crop LAI, and performed more robustly than the reflectance-LUT. In this study, UAV multispectral images were acquired during the wheat grain filling period in different years, which might be an effective way to demonstrate the robustness of LUT method for LAI retrieval. We speculate that because the VIs were expressed as a nonlinear mathematical combination of different bands (Table 7), the differences among sensed objects identified only by G, B, or R reflectance observations would be increased by VI calculations [52]. Previous studies have suggested that VIs can provide information on vegetation phenotyping and reduce some interferences from the soil, atmosphere, and shadows [17,51]. Accordingly, the VI-LUT performed more accurately and robustly than the reflectance-LUT for LAI retrieval. However, the reflectance-LUT strategy is also suitable for LAI retrieval when the spectral bands are not limited [11,18]. Most of the multispectral cameras applied in agriculture, such as the multiSPEC-4C and Micasense RedEdge-M cameras which are commonly used, only include three to six bands; thus, they do not provide sufficient information on spectral reflectance or reduce influences of environmental noises. Therefore, the VI-LUT strategy is a fine alternative method for LAI retrieval and is particularly recommended for precision-agriculture applications with multispectral cameras having limited bands.

4.2. Analyses of LAI Retrieval Performance for Different VI-LUTs

In this study, the differences in LAI retrieval among the four VI-LUTs were not obvious. This is because most of the bands used for the VI calculation were the same, leading to a strong correlation between the VIs. However, the VI-LUTs had specific characteristics regarding LAI estimation. The MCARI2-LUT and NDVI-LUT performed slightly better than the other VI-LUTs in this study (Table 9 and Table 11). We cannot explain the superiority of NDVI and MCARI2 for LAI retrieval based on the results of this study. However, previous studies have indicated that MCARI2 involving the G band is sensitive to vegetation coverage and types, as the G band can indicate canopy chlorophyll content that is strongly related to the LAI [53,54]. Therefore, MCARI2 is a robust index for LAI estimation [48]. Moreover, in order to make the issue—the better performance of NDVI and MCARI2 for LAI retrieval—more convincing, future studies can further investigate the different formulas and the characteristics of different bands for the calculations of VIs.

Moreover, some VIs, such as NDVI, tend to be saturated for high-density vegetation [14,55]. However, as shown in Figure 6, the saturation problem of VIs was not significant in this study. We speculate that the foregoing results were related to the LAI measurements. In this study, the LAI values were relatively low, compared with previous studies in which the LAI was measured using an LAI-2200, LAI-2000, or digital hemispherical photography (DHP) [56,57]. The LAI measured using an LAI-2200, LAI-2000, and DHP included leaves and other green parts such as stems and ears, whereas the LAI measured using an LI-3000C only includes leaves (no other green parts). We compared the LAI values for the LI-3000C and LAI-2200. The results indicated that the LAI measured using the LAI-2200 was higher than using the LI-3000C. The R² value of the regression model was 0.82, indicating that the LAI measurements in this study were feasible. In this study, if the LAI was measured using an LAI-2200, VIs may have underestimated the LAI (with an increased RMSE and MRE); thus, the saturation problem of VIs might occur when the canopy achieves middle-to-high coverage.

4.3. Other Issues Regarding LAI Retrieval Accuracy

Both the VIs and LAI could characterize the functional traits of vegetation canopies, but the relationships between the LAI and VIs were nonlinear and varied between the different vegetation types [58]. The LAI retrieval accuracy for maize was lower than that for wheat for both two-year multispectral observations (please see Appendix A Figure A1 and Appendix A

Table A1. R², RMSE, and MRE for linear regression between simulated LAI and estimated LAI based on reflectance-LUTs/VI-LUTs using multispectral UAV data in 2016 and 2018.

Crops	Year	Reflectance-LUTs				VIs-LUTs
		fRMSE Input	R2	RMSE	MRE	fRMSE Input	R2	RMSE	MRE
Wheat	2016	R,NIR	0.41	2.49	1.23	NDVI	0.70	0.47	0.16
	n = 72	R,NIR,G	0.41	2.49	1.23	MCARI	0.75	0.83	0.34
		R,NIR,E	0.41	2.49	1.23	NRI	0.74	1.04	0.40
	2018	R,NIR	0.42	0.94	0.70	NDVI	0.76	0.44	0.25
	n = 107	R,NIR,B	0.42	0.94	0.70	ARVI	0.74	0.51	0.30
		R,NIR,G	0.42	0.94	0.70	MCARI	0.75	0.38	0.22
		R,NIR, E	0.42	0.94	0.70	NRI	0.75	0.46	0.25
Maize	2016	R,NIR	0.00	1.48	0.42	NDVI	0.68	0.73	0.17
	n = 40	R,NIR,G	0.00	1.48	0.42	MCARI	0.73	0.58	0.16
		R,NIR,E	0.00	1.48	0.42	NRI	0.61	1.56	0.40
	2018	R,NIR	0.00	2.21	0.87	NDVI	0.61	0.62	0.18
	n = 113	R,NIR,B	0.00	2.21	0.87	ARVI	0.59	0.66	0.19
		R,NIR,G	0.00	2.21	0.87	MCARI	0.71	1.19	0.42
		R,NIR, E	0.00	2.21	0.87	NRI	0.62	0.66	0.20

) and one-year hyperspectral observations (please see Appendix A

Table A2. LAI retrieval accuracy (R², RMSE, and MRE) based on reflectance-LUTs/VIs-LUTs using hyperspectral UAV data in 2018 (n = 89).

Crops	Reflectance-LUTs				VIs-LUTs
	Bands	R2	RMSE	MRE	VIs	R2	RMSE	MRE
Wheat	R, NIR	0.27	2.576	2.22	NDVI	0.80	0.55	0.31
	R, NIR, B	0.27	2.576	2.22	ARVI	0.76	0.66	0.40
	R, NIR, G	0.27	2.579	2.22	MCARI2	0.82	0.37	0.26
	R, NIR, E	0.32	2.603	2.24	NRI	0.75	0.58	0.32
Maize	R, NIR	0.07	2.599	2.56	NDVI	0.22	0.75	0.36
	R, NIR, B	0.12	2.550	2.50	ARVI	0.28	0.90	0.45
	R, NIR, G	0.12	2.550	2.50	MCARI2	0.39	0.44	0.17
	R, NIR, E	0.11	2.547	2.50	NRI	0.33	0.78	0.37

). The main reason for this might be the different canopy structures. The remotely sensed LAI implicitly refers to green parts that could be either leaves or other green elements of a vegetation canopy [51]. Thus, the LAI is sensitive to plant structure, such as the clump level, leaf angle, and crop height [59,60]. The PROSAIL model was successfully used for various vegetation types. It treats the canopy as a collection of absorbing and scattering tissues randomly distributed in a horizontal layer. Hence, the PROSAIL model is recommended for application to homogenous crops [11,61]. Both maize and wheat are row crops. When they are at the early stage of growth, they have incomplete coverage and strong leaf clumping, and the background reflectance dominates the spectral signal and affects LAI retrieval [43]. At the middle and later growth stages, their canopies tend to be homogenous. Therefore, the PROSAIL retrieval performance was limited for hyperspectral observations of maize at the early stage of growth (R² ≤ 0.39) (please see Appendix A Table A2), in agreement with previous studies [62]. Furthermore, the biophysical and biochemical characteristics of vegetation change during different stages. In this study, we did not focus on the effects of phenology on the LAI retrieval. However, continuous UAV observations of critical crop growth periods are required for further validations of the robustness and accuracy of VI-LUTs for LAI retrieval.

Statistical regression models (e.g., a simple linear regression model between vegetation parameters and VIs) and machine-learning models (e.g., random forest models and support vector models) have been widely employed for the retrieval of crop biophysical and biochemical parameters. However, one of the main drawbacks is that their applications are valid only in the areas for which they have been calibrated [63]. The PROSAIL model, as one of the radiative transfer models, can overcome this disadvantage and does not require field measurements in conjunction with remote observations for calibration and validation. Because various combinations of canopy parameters may yield similar spectral reflectance, the ill-posed retrieval problem of the PROSAIL model hinders LAI retrieval. Therefore, a priori information, i.e., ranges of model inputs for model parameterization, is essential. The a priori information is related to many factors, such as the crop types, canopy structure, and growth period of the crops. Thus, the low-dimensional VI-LUT should be modified according to the specific vegetation and growth periods.

4.4. Analyses of Optimization for LUT Strategies

In this study, the VI formula structure may also have affected VI calculation in LAI retrieval. Additionally, the central wavelength and bandwidth of VIs affected LAI retrieval. LAI retrieval with four simulated datasets indicate that the VI-LUTs with optimized narrow bands (4 nm) had the highest accuracy. However, the effects of different bandwidths on the accuracy of crop parameter estimation remain unclear. Some hyperspectral narrow bands are sensitive to certain crop parameters, but they might be insensitive to other crop parameters, which may increase the computational loads and distort the accuracy of LAI retrieval [41]. Observational data with broad bands can yield strong reflectance signals, which may reduce the influence of the surrounding environment. Therefore, when selecting the bandwidth in practice, it is essential to achieve a tradeoff between the retrieval sensitivity and the signal-noise ratio.

In summary, this study demonstrates that the VI-LUT strategy based on the PROSAIL model is an alternative method for LAI retrieval with higher accuracy and robustness than the reflectance-LUT strategy. It is valuable for monitoring the crop growth status in precision-agriculture owing to its advantages of simplicity, robustness, accuracy, and saving time and labor. Moreover, the optimization of the central wavelengths and bandwidths of VIs may be helpful for designing multispectral cameras according to specific applications.

5. Conclusions

The objective of this study was to develop a crop LAI retrieval strategy with accuracy, robustness, and ease of use for UAV remote sensing applications in precision-agriculture at the field scale, which would assist the design of UAV multispectral cameras for agronomic monitoring. The study was conducted using UAV remotely sensed images (two-year multispectral data and one-year hyperspectral data) and the PROSAIL model simulation approach. It was concluded that the low-dimensional VI-LUTs were easy to use and more robust (R² ≥ 0.74, RMSE ≤ 0.51, MRE ≤ 0.31 for multispectral datasets) than reflectance-LUTs (R² ≤ 0.42, RMSE ≥ 0.94, MRE ≥ 0.70 for multispectral datasets) in crop LAI retrieval based on UAV data. The differences in LAI retrieval among four VI-LUTs were not obvious. However, the MCARI2-LUT and NDVI-LUT performed slightly better than the other VI-LUTs and thus are recommended for crop LAI retrieval. Moreover, both the central wavelengths and bandwidths of the VIs affected LAI retrieval. The hyperspectral UAV data with central wavelengths of 756 nm for the NIR band and 612 nm for the R band and a narrow bandwidth (~4 nm) improved the performance of the MCARI2-LUT for LAI retrieval.

The VI-LUT strategy based on the PROSAIL model, as an alternative method for LAI retrieval at the field scale, is recommended for crop growth monitoring with multispectral cameras having limited bands. The optimized central wavelengths and bandwidths of VIs, as well as corresponding methods of VI optimization, might contribute to improving the design of multispectral cameras onboard UAVs for retrieving vegetation traits in precision-agriculture.