An Approach for Joint Estimation of Grassland Leaf Area Index and Leaf Chlorophyll Content from UAV Hyperspectral Data

Abstract

Leaf area index (LAI) and leaf chlorophyll content (Cab) are two important indicators of vegetation growth. Due to the high-coupling of spectral signals of leaf area and chlorophyll content, simultaneous retrieval of LAI and Cab from remotely sensed date is always challenging. In this paper, an approach for joint estimation of grassland LAI and Cab from unmanned aerial vehicle (UAV) hyperspectral data was proposed. Firstly, based on a PROSAIL model, 15 typical hyperspectral vegetation indices (VIs) were calculated and analyzed to identify optimal VIs for LAI and Cab estimation. Secondly, four pairs of VIs were established and their discreteness was also calculated for building a two-dimension matrix. Thirdly, a two-layer VI matrix was generated to determine the relationship of VIs with LAI values and Cab values. Finally, LAI and Cab were jointly retrieved according to the cells of the two-layer matrix. The retrieval reduced the cross-influence between LAI and Cab. Compared with the VI empirical model and the single-layer VI matrix, the accuracy of LAI and Cab retrieved from UAV hyperspectral data based on the two-layer VI matrix was significantly improved (for LAI: R² = 0.73, RMSE = 0.91 m²/m² and u(SD) = 0.82 m²/m²; for Cab: R² = 0.79, RMSE = 11.7 μg/cm² and u(SD) = 10.84 μg/cm²). The proposed method has the potential for rapid retrieval of LAI and Cab from hyperspectral data. As a method similar to look-up table, the two-layer matrix can be used directly for LAI and Cab estimation without the need for prior measurements for training.

1. Introduction

Grassland degradation is a serious ecological and economical problem in northern grassland of China. Leaf area index (LAI) and leaf chlorophyll content (Cab), two key indicators of grassland vegetation growth, are important access points for studying the problem [1,2,3]. Due to the mobility characteristics of unmanned aerial vehicles (UAVs), they provide a more efficient and convenient means for grasslands monitoring [4,5,6,7]. Relatively high spatio-temporal resolution and spectral resolution information for grasslands can be collected based on the UAV system, which can be equipped with consumer-grade hyperspectral sensors.

Modeling and analysis based on hyperspectral data are important methods to estimate LAI and Cab in large areas of grassland. Many studies have successfully retrieved LAI and Cab using the inversion of a radiative transfer model using methods such as look-up tables (LUT), neural networks (NNs), simulated annealing (SA), genetic algorithms (GAs) and other optimization techniques [2,8,9,10,11]. However, due to the measurement and model uncertainties and the limited information carried by the radiometric signal, the inversion of physical models is ill-posed [12,13,14]. Vegetation indices (VIs), calculated based on target information of spectral signals, are computationally simple. VIs can be used as the model parameter in place of the full spectral band for inversion model construction, representing an important method for retrieving LAI and Cab from remote sensing data.

The changes in LAI and Cab are jointly reflected in canopy reflectance. At times, two different vegetation types or growth states (different combinations of LAI and Cab) may present the same spectral characteristics. Therefore, under the coupling effect of spectral information, the retrieval of LAI and Cab based on VIs is influenced by their variation [15]. VI optimization and VI combinations are two efficient ways to improve the retrieval accuracy of LAI and Cab. The optimized VIs, such as TCARI and MCARI, have been adjusted to improve their relevance to Cab by reducing the influence of background soil spectral information [16,17]. However, the effects of LAI on Cab retrieval still exist. Liang et al. pointed out that although the OSAVI had become more sensitive to LAI through adjustment, the Cab still had influence on the relationships between the VI and LAI at low chlorophyll contents [15]. VI combinations have been established to reduce the effect of one parameter on the retrieval of another parameter. For example, since VIs such as TCARI and OSAVI are sensitive to Cab and LAI, respectively, their ratio TCARI/OSAVI can be used for Cab retrieval by reducing the effect of LAI [16]. The simple ratio indices include R750/R700 and R865/R665; the former is strongly correlated to Cab, while the latter is widely used in LAI estimation. Therefore, their combination can be used for the estimation of Cab and LAI [17].

Simultaneous retrieval of LAI and Cab based on appropriate VI combinations is one way to solve the problem and improve the estimation accuracy. In this study, an approach for the joint estimation of grassland LAI and Cab from unmanned aerial vehicle (UAV) hyperspectral data is proposed. The approach focuses on (1) evaluating different VI combinations for joint estimation of LAI and Cab based on PROSAIL model simulation; (2) generating a two-layer VI matrix to determine the relationship between VIs and LAI values and Cab values; (3) testing the sensitivities and reliability of the two-layer matrix with simulated data; (4) carrying out the retrieval of LAI and Cab from UAV hyperspectral data over a grassland demonstration area.

2. Materials and Methods

2.1. The Study Area and Data

2.1.1. Description of the Study Area and In Situ Measurements

In order to carry out research on grassland ecosystem function assessment, a comprehensive field campaign was conducted over Inner Mongolia Grassland Ecosystem Research Station (IMGERS). The IMGERS (UL: 43°34′18.81″N, 116°39′24.82″E; UR: 43°34′19.03″N, 116°40′40.75″E; BL: 43°33′7.34″N, 116°39′25.23″E; BR: 43°33′7.56″N, 116°40′41.13″E) is located on the north side of the Xilin river basin, which is one of the most representative steppe zones in China, as shown in Figure 1. The Xilin river basin is located between 43°26′–44°29′N and 115°32′–117°12′E in the Inner Mongolia Autonomous Region, China. In the Xilin river basin, the annual mean temperature is about 1–4 °C, while annual mean precipitation is ~350 mm, most of which is concentrated between June and August. Stipa grandis and Leymus chinensis are two dominant species in the basin [18,19].

During the campaign, surface reflectance of soil, grassland leaf and grassland canopy was collected on 10 Aug 2018 using SVC HR-1024 manufactured by Spectra Vista Corporation, Poughkeepsie, NY, USA. The SVC HR-1024 is a field portable spectrometer, and its spectral range covers from 0.35 μm to 2.5 μm. During the field campaign, surface reflectance of each ESU (Elementary Sampling Unit, 1 m × 1 m) was calculated using the average of 5 measurements.

A Plant Canopy Analyzer LAI-2200 was used for collecting in situ LAI shortly after sunrise or shortly before sunset. The average LAI of each ESU was calculated based on four low-canopy measurements and one above-canopy measurement. After LAI measurements, top, middle and bottom leaves of the plant were cut and used for chlorophyll content (Cab) measurement in the laboratory.

2.1.2. Description of the UAV Hyperspectral Data

A 720 g Rikola hyperspectral camera was carried on the UAV platform and used for obtaining the hyperspectral data over the IMGERS. The Rikola hyperspectral camera was a full-frame lightweight device, provided by SenopOy, Kangasala, Finland.

Table 1. Technical properties of the Rikola hyperspectral camera.

Parameter	Value
Image mode	Frame-based
FOV	36.5°
Pixel resolution	1010 × 1010
Field resolution	6.5 cm@100 m
Spectral range	500–900 nm
Spectral resolution	~10 nm
Weight	720 g

lists the technical properties of the Rikola hyperspectral camera [20].

The HSIs required a series of pre-processing steps, including band matching, geometric, radiometric and atmospheric corrections. The band matching and full image matching were processed using supporting software RegMosaic. Georeferencing of HSIs scenes to the reference system WGS 84 UTM 50 N was carried out based on 13 field geometric control points. Since the atmospheric influence was almost negligible, an empirical line approach was carried out for the atmospheric correction and the conversion from radiance to reflectance of the HSI spectra. During the processing, the PVC panels in white, grey and black were taken as the spectrally defined ground targets for calibration.

2.2. A Two-Layer Matrix Approach for Joint Estimation of LAI and Cab

2.2.1. PROSAIL Model and Simulated Dataset

The PROSAIL model, coupling PROSPECT leaf optical properties model [21] and SAIL canopy bidirectional reflectance model [22], has been used for more than twenty years to study plant canopy spectral and directional reflectance in the solar domain and vegetation biophysical properties [23,24]. In this article, the PROSAIL model was used for the selection and evaluation of vegetation indices. The PROSPECT model simulates the leaf’s optical properties (the reflectance and transmittance) from 400 nm to 2500 nm as a function of leaf structure parameters and biochemical parameters [24], including equivalent water thickness C_w (g/cm²), dry matter content C_m (g/cm²), leaf chlorophyll a and b content C_ab (μ/cm²) and leaf structure parameter N (unitless). The reflectance and transmittance of leaves simulated by PROSPECT are then inputted into SAIL model to simulate the canopy spectral and directional reflectance as a function of canopy structure parameters [24], including soil reflectance ρ_s (unitless), leaf area index LAI (m²/m²), average leaf angle ALA (°), hot-spot size parameter Hot (m/m) and external parameters for view zenith angle VZA (°), sun zenith angle SZA (°) and relative azimuth angle RSA (°) between the sensor and sun.

An appropriate parameterization of the PROSAIL model was required for simulating canopy reflectance. Based on the field measurements and documents [23,24,25,26,27,28], the ranges of input parameters are shown in Table 2. Background reflectance can be a significant contributor to the canopy reflectance signal [28]. In this study, a factor f was established based on field measured soil spectrum for quantizing the variation in soil brightness.

$$soil=f\times soi{l}_{min}+\left(1-f\right)\times soi{l}_{max}$$

(1)

where $soi{l}_{min}$ and $soi{l}_{max}$ are the minimum and maximum of field measured soil reflectance spectra.

Table 2. The main input parameters of the PROSAIL model.

Parameters	Symbol	Unit	Expectation	Range
Chlorophyll a + b content	Cab	μg/cm2	45	10–90
Dry matter content	Cm	g/cm2	0.005	0.001–0.01
Equivalent water thickness	Cw	g/cm2	0.02	0.005–0.04
Leaf area index	LAI	m2/m2	2.8	0–5
Average leaf angle	ALA	°	15	10–65
Leaf structure parameter	N		1.5	1.4–2.2
Hot-spot size parameter	Hot	m/m	0.2	0–1.4
Soil factor	f		0.4	0.1–0.9
View zenith angle	VZA	°		5.78
Sun zenith angle	SZA	°		23.12
Relative azimuth angle	RSA	°		111.39

Table 2. The main input parameters of the PROSAIL model.

Parameters	Symbol	Unit	Expectation	Range
Chlorophyll a + b content	Cab	μg/cm2	45	10–90
Dry matter content	Cm	g/cm2	0.005	0.001–0.01
Equivalent water thickness	Cw	g/cm2	0.02	0.005–0.04
Leaf area index	LAI	m2/m2	2.8	0–5
Average leaf angle	ALA	°	15	10–65
Leaf structure parameter	N		1.5	1.4–2.2
Hot-spot size parameter	Hot	m/m	0.2	0–1.4
Soil factor	f		0.4	0.1–0.9
View zenith angle	VZA	°		5.78
Sun zenith angle	SZA	°		23.12
Relative azimuth angle	RSA	°		111.39

With the defined ranges, 150,000 parameter combinations were obtained using uniform random sampling to generate 150,000 canopy spectra. The simulated dataset was divided into two parts: 100,000 simulations were randomly selected and used for modeling and the other simulations were used independently for model testing.

2.2.2. Selection of the Vegetation Indices

Based on the simulated canopy spectral data, 15 typical vegetation indices were calculated and analyzed for LAI and Cab estimation, including normalized difference vegetation index (NDVI), simple ratio vegetation index (SR), etc. The 15 typical vegetation indices are listed in Table 3. The wavelengths used were within the HSI’s band domains, in the visible and near-infrared ranges.

Table 3. The 15 typical vegetation indices.

Vegetation Index	Formulation	Reference
NDVI	$\left({\rho }_{864}-{\rho }_{664}\right)/\left({\rho }_{864}+{\rho }_{664}\right)$	[29]
SR	${\rho }_{752}/{\rho }_{704}$	[30]
EVI	$2.5\times \frac{{\rho }_{864}-{\rho }_{664}}{{\rho }_{864}+6\times {\rho }_{664}-7.5\times {\rho }_{504}+1}$	[31]
MTVI1	$1.2\left[1.2\left({\rho }_{800}-{\rho }_{550}\right)-2.5\left({\rho }_{672}-{\rho }_{550}\right)\right]$	[28]
MTVI2	$\frac{1.5\left[1.2\left({\rho }_{800}-{\rho }_{550}\right)-2.5\left({\rho }_{672}-{\rho }_{550}\right)\right]}{\sqrt{{\left(2{\rho }_{800}+1\right)}^2-\left(6{\rho }_{800}-5\sqrt{{\rho }_{672}}\right)-0.5}}$	[28]
OSAVI	$\frac{\left(1+0.16\right)\times \left({\rho }_{800}-{\rho }_{672}\right)}{\left({\rho }_{800}+{\rho }_{672}+0.16\right)}$	[32]
MSAVI	$\frac{2{\rho }_{800}+1-\sqrt{(2{\rho }_{800}+1{)}^2-8\left({\rho }_{800}-{\rho }_{672}\right)}}{2}$	[33]
TCARI	$3[({\mathsf{\rho }}_{704}-{\mathsf{\rho }}_{672})-0.2({\mathsf{\rho }}_{704}-{\mathsf{\rho }}_{552})\left(\frac{{\mathsf{\rho }}_{704}}{{\mathsf{\rho }}_{672}}\right)]$	[17]
TCARI2	$3\left[({\mathsf{\rho }}_{752}-{\mathsf{\rho }}_{704})-0.2({\mathsf{\rho }}_{752}-{\mathsf{\rho }}_{552}){\mathsf{\rho }}_{752}/{\mathsf{\rho }}_{704}\right]$	[34]
MCARI	$\left[{\mathsf{\rho }}_{704}-{\mathsf{\rho }}_{672}-0.2\left({\mathsf{\rho }}_{704}-{\mathsf{\rho }}_{552}\right)\right]\times \left({\mathsf{\rho }}_{704}/{\mathsf{\rho }}_{672}\right)$	[35]
MTCI	$\left({\mathsf{\rho }}_{752}-{\mathsf{\rho }}_{712}\right)/\left({\mathsf{\rho }}_{712}-{\mathsf{\rho }}_{680}\right)$	[36]
CIgreen	${\rho }_{864}/{\rho }_{560}$−1	[37]
CIrededge	${\rho }_{864}/{\rho }_{720}$−1	[37]
REP	$704+35\frac{\left({\rho }_{672}+{\rho }_{784}\right)/2-{\mathsf{\rho }}_{704}}{{\rho }_{744}-{\mathsf{\rho }}_{704}}$	[38]
TVI	$0.5\left[120\left({\rho }_{752}-{\rho }_{552}\right)-200\left({\rho }_{672}-{\rho }_{552}\right)\right]$	[39]

Table 3. The 15 typical vegetation indices.

Vegetation Index	Formulation	Reference
NDVI	$\left({\rho }_{864}-{\rho }_{664}\right)/\left({\rho }_{864}+{\rho }_{664}\right)$	[29]
SR	${\rho }_{752}/{\rho }_{704}$	[30]
EVI	$2.5\times \frac{{\rho }_{864}-{\rho }_{664}}{{\rho }_{864}+6\times {\rho }_{664}-7.5\times {\rho }_{504}+1}$	[31]
MTVI1	$1.2\left[1.2\left({\rho }_{800}-{\rho }_{550}\right)-2.5\left({\rho }_{672}-{\rho }_{550}\right)\right]$	[28]
MTVI2	$\frac{1.5\left[1.2\left({\rho }_{800}-{\rho }_{550}\right)-2.5\left({\rho }_{672}-{\rho }_{550}\right)\right]}{\sqrt{{\left(2{\rho }_{800}+1\right)}^2-\left(6{\rho }_{800}-5\sqrt{{\rho }_{672}}\right)-0.5}}$	[28]
OSAVI	$\frac{\left(1+0.16\right)\times \left({\rho }_{800}-{\rho }_{672}\right)}{\left({\rho }_{800}+{\rho }_{672}+0.16\right)}$	[32]
MSAVI	$\frac{2{\rho }_{800}+1-\sqrt{(2{\rho }_{800}+1{)}^2-8\left({\rho }_{800}-{\rho }_{672}\right)}}{2}$	[33]
TCARI	$3[({\mathsf{\rho }}_{704}-{\mathsf{\rho }}_{672})-0.2({\mathsf{\rho }}_{704}-{\mathsf{\rho }}_{552})\left(\frac{{\mathsf{\rho }}_{704}}{{\mathsf{\rho }}_{672}}\right)]$	[17]
TCARI2	$3\left[({\mathsf{\rho }}_{752}-{\mathsf{\rho }}_{704})-0.2({\mathsf{\rho }}_{752}-{\mathsf{\rho }}_{552}){\mathsf{\rho }}_{752}/{\mathsf{\rho }}_{704}\right]$	[34]
MCARI	$\left[{\mathsf{\rho }}_{704}-{\mathsf{\rho }}_{672}-0.2\left({\mathsf{\rho }}_{704}-{\mathsf{\rho }}_{552}\right)\right]\times \left({\mathsf{\rho }}_{704}/{\mathsf{\rho }}_{672}\right)$	[35]
MTCI	$\left({\mathsf{\rho }}_{752}-{\mathsf{\rho }}_{712}\right)/\left({\mathsf{\rho }}_{712}-{\mathsf{\rho }}_{680}\right)$	[36]
CIgreen	${\rho }_{864}/{\rho }_{560}$−1	[37]
CIrededge	${\rho }_{864}/{\rho }_{720}$−1	[37]
REP	$704+35\frac{\left({\rho }_{672}+{\rho }_{784}\right)/2-{\mathsf{\rho }}_{704}}{{\rho }_{744}-{\mathsf{\rho }}_{704}}$	[38]
TVI	$0.5\left[120\left({\rho }_{752}-{\rho }_{552}\right)-200\left({\rho }_{672}-{\rho }_{552}\right)\right]$	[39]

To identify appropriate vegetation indices for LAI and Cab estimation, the VIs shown in Table 2 were tested and analyzed based on the performance of different fitting models. Firstly, the VIs were taken as the independent variable, vegetation parameter (LAI or Cab) was taken as the dependent variable. Then, the optimal fitting models were constructed based on the best fit from the linear regression, power regression, logarithmic regression and exponential regression. Finally, the coefficient of determination (R²) and root mean square error (RMSE) of the curve fitting models were used for assessing the performance of VIs.

The inversion performance of the 15 VIs for LAI and Cab was assessed and the RMSE/R² values are ranked in descending order in

Table 4. The relationship between LAI and VIs.

Indices	Fitting Equation	a	b	R2	RMSE
NDVI	f(x) = axb	3.133	2.036	0.7442	0.4395
OSAVI	f(x) = axb	3.633	1.544	0.7325	0.4494
MTVI2	f(x) = axb	3.06	0.87	0.7081	0.4695
MSAVI	f(x) = ax + b	3.574	−0.1135	0.683	0.4892
EVI	f(x) = ax + b	3.421	−0.2812	0.6279	0.5301
MTVI1	f(x) = axb	3.02	0.8121	0.6277	0.5302
TVI	f(x) = axb	0.1607	0.8104	0.6225	0.5339
SR	f(x) = axb	0.6915	0.7908	0.4564	0.6407
CIgreen	f(x) = axb	0.8082	0.5259	0.4449	0.6474
MCARI	f(x) = axb	3.035	0.3223	0.4306	0.6557
CIrededge	f(x) = axb	1.937	0.599	0.3112	0.7212
TCARI	f(x) = ax + b	4.22	0.9238	0.1814	0.7862
MTCI	f(x) = aexp(bx)	1.025	0.2256	0.0954	0.8265
TCARI2	f(x) = aexp(bx)	1.592	−0.4836	0.0111	0.8641
REP	f(x) = aexp(bx)	0.0018	0.0093	0.0061	0.8663

and

Table 5. The relationship between Cab and VIs.

Indices	Fitting Equation	a	b	R2	RMSE
REP	f(x) = axb	2.142 × 10−210	73.98	0.8104	8.658
MTCI	f(x) = axb	31.85	0.7857	0.7199	11.133
CIrededge	f(x) = axb	58.03	0.5938	0.5107	14.23
TCARI	f(x) = ax + b	−159.1	66.56	0.4704	14.81
CIgreen	f(x) = ax + b	4.722	27.67	0.3345	16.6
MCARI	f(x) = aexp(bx)	59.5	−2.017	0.3265	16.7
SR	f(x) = ax + b	8.852	20.57	0.3018	17
TCARI2	f(x) = ax + b	−71.88	54.47	0.2559	17.55
NDVI	f(x) = aexp(bx)	35.4	0.3531	0.0208	20.13
TVI	f(x) = aexp(bx)	48.94	−0.0053	0.0097	20.25
MTVI1	f(x) = ax + b	−8.827	48.79	0.0093	20.25
EVI	f(x) = aexp(bx)	40.37	0.2029	0.0078	20.27
OSAVI	f(x) = aexp(bx)	40.06	0.2076	0.0072	20.27
MSAVI	f(x) = aexp(bx)	41.33	0.1824	0.0062	20.28
MTVI2	f(x) = axb	43.87	−0.0226	0.0018	20.33

.

Based on the above results, NDVI (R² = 0.7442) and OSAVI (R² = 0.7325) were selected for LAI estimation, while REP (R² = 0.8104) and MTCI (R² = 0.7199) were selected for Cab estimation. The different combinations of these vegetation indices were established as shown in

Table 6. The VI combinations.

↓ LAI	REP	MTCI
NDVI	(NDVI, REP)	(NDVI, MTCI)
OSAVI	(OSAVI, REP)	(OSAVI, MTCI)
↓ Cab	NDVI	OSAVI
REP	(REP, NDVI)	(REP, OSAVI)
MTCI	(MTCI, NDVI)	(MTCI, OSAVI)

.

2.2.3. Evaluation of VI Combinations

An appropriate VI combination for vegetation parameter estimation should be not only sensitive to the target parameters, but also insensitive to interference factors [40]. To evaluate the sensitivity of the vegetation index combinations to different vegetation parameters, the PROSAIL model was used to simulate the influence of LAI, Cab, C_m, C_w, N, ALA and ρ_s on the 15 vegetation indices described above. During the sensitivity analysis process, only the interference factor was varied according to the threshold range and the remained parameters of PROSAIL model were set to the fixed value. The VI combinations were drawn in a two-VI space, as shown in Figure 2.

In Figure 2, the LAI isolines (the same symbol) and the Cab isolines (the same color) had different trends. This indicates that the LAI and Cab can be well-separated by VI combinations. In general, the sensitivity of VI combinations changed with the LAI and Cab, varying in different VI spaces and positions. The VI combinations of NDVI–MTCI and OSAVI–MTCI had higher discreteness (the Cab isolines in these VI spaces were the most separated) and are better for distinguishing Cab, indicating the most sensitivity. Due to the saturation effects, the higher the LAI value, the lower the sensitivity of the VI combinations to LAI and the more dense the LAI isolines in two-VI space. When both the LAI and Cab were high, the saturation effect was most pronounced.

It was difficult to directly judge the dispersion of VI combination only from the figure. Therefore, the ability of different combinations to distinguish LAI or Cab was quantitatively described based on the distance of adjacent points in the matrix. The operations were as follows: (1) Normalize the matrix; (2) Calculate the mean value (Lave) and standard deviation (Lstd) of the distance between different Cab points on the LAI isoline to represent the dispersion of the current matrix to Cab. The larger the Lave, the better the dispersion of the matrix to Cab, and the smaller the Lstd, the better the stability of the matrix. Similarly, the distance mean and variance of LAI on Cab isolines were also calculated. The larger the Lave of LAI, the better the dispersion of the matrix for LAI, and the smaller the Lstd of LAI, the better the stability of the matrix.

As shown in

Table 7. The degree of dispersion of different matrices.

	NDVI–REP	OSAVI–REP	NDVI–MTCI	OSAVI–MTCI
Lave of Cab	0.0886	0.0849	0.0863	0.0931
Lstd of Cab	0.0946	0.0888	0.0359	0.0285
Lave of LAI	0.0694	0.0731	0.0716	0.0683
Lstd of LAI	0.0971	0.0819	0.0896	0.0897

, for Cab, OSAVI–MTCI not only had the highest dispersion (Lave of Cab was 0.0931) but also the lowest standard deviation (Lstd of Cab was 0.0285), which was consistent with the phenomenon described above. For LAI, OSAVI–REP performed best, but at the same time, it had certain instability for Cab. It can also be seen from Figure 2 that in the OSAVI–REP space, contour lines such as Cab = 5 μg/cm², 20 μg/cm², 35 μg/cm² were relatively discrete, while contour lines such as Cab = 65 μg/cm², 80 μg/cm², 90 μg/cm² were relatively clustered. This indicates that the degree of dispersion variability of this VI combination is relatively large (Lstd of Cab was 0.0888). Thus, a partitioned multi-stage inversion strategy should be established for joint estimation of LAI and Cab. For example, the Cab parameter should be preferentially retrieved because its VI combinations were more sensitive to Cab.

The dispersion of OSAVI–REP to Cab under different thresholds was further analyzed to determine the combination mode of the two VI matrices. It can be seen from Figure 3 that the dispersion degree of OSAVI–REP decreases significantly with the increase in Cab value. Therefore, it is necessary to introduce the VI matrix in stages according to the variation trend of the dispersion. That is, OSAVI–REP should be used when the Cab is low, while OSAVI–MTCI should be used when the Cab is high.

The effects of the other parameters on VI combinations are also evaluated in this paper. Taking VI combination OSAVI–MTCI as an example, its sensitivity to the parameters of ALA, Cm, Cw, N and ρ_s was analyzed using the simulations of PROSAIL. The results are shown in Figure 4. By comparing the curves in Figure 4 with those of the curves in Figure 2, the variation in Cm, Cw and ρ_s had no obvious effect on the results. The changes in ALA and N had limited effects on OSAVI–MTCI combination, for which N had a relatively high impact; the parameter ALA had a great influence on LAI when LAI value was low (Figure 4(a1)–(a4)), but it did little to change the overall properties of OSAVI–MTCI combination. Although the effect of N on Cab was minimal when the Cab value was lower than 35 μg/cm² (Figure 4(d1)–(d4)), it led to a reduction in the dispersion of OSAVI–MTCI combination (the Cab curves were relatively gathered).

2.2.4. Establishment of a Two-Layer Matrix Inversion Method

As shown in Figure 5, the VI combinations of OSAVI–REP and OSAVI–MTCI were divided into cells to establish the inversion matrix. Each matrix had 1000 × 1000 cells and each cell corresponded to a small range of vegetation index values. Based on the modeling dataset simulated from the PROSAIL model, the LAI and Cab values were assigned to a matrix cell according to the values of the two simulated VIs. Because one cell had multiple values of LAI and Cab, all the values mapped to the same cell would be calculated as average and standard deviation. Then, the average (ave) and standard deviation (sd) values of the LAI and Cab were mapped to each cell, constituting the average and standard deviation matrices.

The distribution of matrix dispersion was uneven due to the value change of LAI and Cab, which was described in Section 2.2.3. The dispersion curves of OSAVI–REP matrix and OSAVI–MTCI matrix converge at Cab ≈ 50 μg/cm² (Figure 3). When Cab was low, OSAVI–REP matrix had better discretization, which should be used for Cab estimation. Otherwise, OSAVI–MTCI should be used for Cab estimation. In order to reduce the error brought by the instability of the dispersion of the matrix, the range (40 μg/cm², 60 μg/cm²) was determined as the intermediate choice. Then, the weight coefficient can be established to determine the participation of matrix cells in different layers. As described above, when Cab value ≤ 40 μg/cm², OSAVI–REP dispersion was high. Therefore, the corresponding cell in OSAVI–REP matrix should be assigned a higher weight, which was set to 1 in this paper. Similarly, when Cab value > 60 μg/cm², OSAVI–MTCI dispersion was high and the corresponding cell in OSAVI–MTCI should be assigned a higher weight. When 40 μg/cm² < Cab value ≤ 60 μg/cm², the corresponding cell in OSAVI–REP matrix and OSAVI–MTCI matrix were both assigned as 0.5. At last, the weight coefficient (w) was entered in each cell to form the weight matrix.

In the retrieval process, the expected estimation value of LAI and Cab were calculated based on the average matrix and the weight matrix. Additionally, the total standard deviation represented the inversion uncertainty.

After establishing the two-layer matrix, the LAI and Cab can be easily retrieved based on the four VIs calculated from the reflectance data. If the cell value of each layer matrix was valid, its value can be calculated with the weight coefficient and set as the estimated LAI and Cab. However, in some cases, the VIs calculated from the remotely sensed data were not within the range of the matrix and LAI and Cab values of the corresponding units could not be searched. In other words, the VI combinations calculated based on simulated reflectance data still cannot cover those calculated based on the remotely sensed data. This could be due to improper parameter settings during PROSAIL model simulation or the observation deviation of remotely sensed data. If the searched LAI and Cab values in the located cell were both null in two-layer matrices, the nearest 8-connected neighborhood was averaged and used as the estimated LAI and Cab. If the values of LAI and Cab in the located cell were null only in a one-layer matrix, the non-null values of the other layer were directly set to the estimated LAI and Cab.

2.2.5. Joint Estimation of LAI and Cab

The joint estimation of LAI and Cab included two cases:

Case 1 based on simulated data: the simulated test datasets were used for retrieving LAI and Cab through VI empirical model, the single-layer VI matrix and the two-layer matrix.

As described in Section 2.2, 50,000 canopy spectra were used as a test dataset (TD data) for joint estimation of LAI and Cab. In this case, to further consider bias in the acquisition and processing of remote sensing data, a new test dataset (NTD data) was constructed by adding 5% relative Gaussian noise to the canopy spectra as reflectance uncertainty. The new test dataset was produced, which can also be used for evaluating the influence of reflectance measurement uncertainty on the retrieval accuracy of vegetation parameters.

Therefore, in this case, the two-layer matrix was evaluated by retrieving LAI and Cab from the simulated data with and without noise.

Case 2 based on UAV hyperspectral data: the UAV hyperspectral data obtained during the field campaign were used for retrieving LAI and Cab through approaches of matrices and machine learning, including the single-layer matrix of OSAVI–REP, the single-layer matrix of OSAVI–MTCI, the two-layer matrix, partial least squares regression (PLSR) and random forest (RF).

In Case 2, to improve the model robustness when the model was used for actual observation data, the two-layer matrix was optimized by adding 5% relative Gaussian noise for comparison. Therefore, in this case, the two-layer matrix established with and without noise was evaluated by retrieving LAI and Cab from the UAV hyperspectral data.

3. Results

3.1. Evaluation of the Two-Layer Matrix Based on Simulated Data

Based on the two-layer matrix, LAI and Cab were jointly estimated from the simulated TD and NTD datasets. The inversion results were compared with those estimated using the VI empirical model and the single-layer VI matrix (

Table 8. The retrieving results from simulated data.

Retrieval Methods	Performance (u(SD))		TD Data		NTD Data
			R2	RMSE	R2	RMSE
relationships of LAI–OSAVI: y = 3.633x1.544	LAI	-	0.69	1.06	0.58	1.22
relationships of Cab–REP: y = 2.142 × 10−210x73.98	Cab	-	0.75	12.63	0.69	14.02
single-layer matrix of OSAVI–REP	LAI	0.87	0.75	0.94	0.72	0.99
	Cab	9.51	0.79	12.24	0.73	13.41
single-layer matrix of OSAVI–MTCI	LAI	0.98	0.71	1.01	0.66	1.09
	Cab	7.93	0.83	14.52	0.77	16.16
two-layer matrix of OSAVI–REP and OSAVI–MTCI	LAI	0.67	0.79	0.87	0.76	0.92
	Cab	7.6	0.85	11.05	0.79	12.66

). The following statistical indicators were used for accuracy analysis: coefficient of determination (R²) and root mean square error (RMSE). The standard deviation (u(SD)) characterized the performance of the matrices and was also used to represent the uncertainty of LAI and Cab estimation.

As shown in Table 8, regardless of whether the test dataset contains noise, the two-layer retrieval matrix achieved the best estimation accuracies of the LAI and Cab. LAI retrieved from NTD (TD) data obtained R² of 0.76 (0.79) and RMSE of 0.92 m²/m² (0.87 m²/m²). Cab retrieved from NTD (TD) data obtained R² of 0.79 (0.85) and RMSE of 12.66 μg/cm² (11.05 μg/cm²). Although the single-layer matrix of OSAVI–REP had similar estimated accuracy for retrieving LAI from the simulated reflectance with and without reflectance noise, for Cab estimation, the two-layer matrix of OSAVI–REP and OSAVI–MTCI performed better. In addition, compared with the single-layer matrix of OSAVI–REP, the uncertainty of the two-layer matrix applied to Cab estimation was lower. The reason why the two-layer matrix was superior to the single-layer matrix was that the two-layer matrix was combined with the most separable part of the single-layer matrix, which reduced the cross-influence of LAI and Cab and improved the robustness of the inversion process.

Compared with the empirical model, the matrix-based approaches were obviously more accurate for estimating LAI and Cab. By combining OSAVI–REP and OSAVI–MTCI matrices, the two-layer matrix achieved a higher estimation than one-dimensional relationships. Retrieval results based on matrix models were better because the matrices better characterize the relationship between reflectance and LAI and Cab. Empirical models based on VI cannot represent the complex relationship between reflectance and LAI and Cab well.

The accuracy of LAI and Cab retrieved from reflectance data containing noise were significantly lower than that from the reflectance data without noise. The result indicates that the uncertainty of remotely sensed observations had an important effect on the estimation accuracy of vegetation parameters.

3.2. Evaluation of the Two-Layer Matrix Based on UAV Hyperspectral Data

The retrieval methods were applied to the UAV hyperspectral data for retrieving the grassland LAI and Cab. Figure 6 shows the results retrieved from UAV hyperspectral data based on different inversion matrices. The statistical parameters of R², RMSE and u(SD) were used for estimation accuracy and uncertainty analysis.

Among the three inversion matrices of single-layer matrix of OSAVI–REP, single-layer matrix of OSAVI–MTCI and two-layer matrix of OSAVI–REP and OSAVI–MTCI, the estimation using the two-layer matrix yielded the highest accuracy and lowest uncertainty. The R², RMSE and u(SD) for LAI were 0.67, 1.03 m²/m² and 0.9 m²/m² (Figure 6c top). The R², RMSE and u(SD) for Cab were 0.70, 14.2 μg/cm² and 12.36 μg/cm² (Figure 6c bottom). The matrix of OSAVI–MTCI had the lowest accuracy for LAI estimation, with R² = 0.55 and RMSE = 1.26 m²/m² (Figure 6a top). The matrix of OSAVI–REP had the lowest accuracy for Cab, with R² = 0.52 and RMSE = 17.26 μg/cm² (Figure 6a top).

As described in Section 3.1, retrieval accuracy is influenced by the uncertainty of canopy reflectance. Therefore, the matrices were optimized by adding noise for comparison. The performance of the matrices with reflectance noise was also analyzed based on the values of mean standard deviation (u(SD)), which was also used to represent the uncertainty of LAI and Cab estimation when using those approaches.

Table 9. The uncertainty of the matrices.

Retrieval Methods	Performance (u(SD))
single-layer matrix of OSAVI–REP	LAI: 0.88, Cab: 9.34
single-layer matrix of OSAVI–MTCI	LAI: 1.04, Cab: 7.83
two-layer matrix of OSAVI–REP and OSAVI–MTCI	LAI: 0.72, Cab: 7.4

shows that the uncertainty of the matrix did not increase significantly after the addition of simulated noise.

The optimized two-layer matrix was applied to UAV hyperspectral data for retrieving the grassland LAI and Cab. The retrieval results were compared with the single-layer matrix of OSAVI–REP and the single-layer matrix of OSAVI–MTCI, as shown in

Table 10. The retrieval of LAI and Cab from UAV hyperspectral data using different methods.

Retrieval Methods	Retrieval Accuracy
	R2	RMSE
single-layer matrix of OSAVI–REP	LAI: 0.66	LAI: 1.04 m2/m2
	Cab: 0.61	Cab: 16.82 μg/cm2
single-layer matrix of OSAVI–MTCI	LAI: 0.63	LAI: 1.13 m2/m2
	Cab: 0.69	Cab: 17.43 μg/cm2
two-layer matrix of OSAVI–REP and OSAVI–MTCI	LAI: 0.73	LAI: 0.91 m2/m2
	Cab: 0.79	Cab: 11.70 μg/cm2
partial least squares regression (PLSR)	LAI: 0.65	LAI: 1.07 m2/m2
	Cab: 0.66	Cab: 14.97 μg/cm2
random forest (RF)	LAI: 0.68	LAI: 1.06 m2/m2
	Cab: 0.70	Cab: 14.88 μg/cm2

. In order to evaluate the applicability of the two-layer matrix more comprehensively, partial least squares regression (PLSR) and random forest (RF) were also selected for comparison. PLSR integrates the ideas of principal component analysis, correlation analysis and multiple regression, and has great advantages in dealing with problems where the number of independent variables is much larger than the number of samples. The RF algorithm is a versatile machine learning method that has a better fitting effect for nonlinear data.

For the three inversion matrices established based on the reflectance dataset with noise, the estimation using the two-layer matrix yielded the highest accuracy and lowest uncertainty. The R², RMSE and u(SD) for LAI were 0.73, 0.91 m²/m² and 0.82 m²/m². The R², RMSE and u(SD) for Cab were 0.79, 11.7 μg/cm² and 10.84 μg/cm². The retrieval results of the matrices with reflectance noise were obviously better than those of matrices without noise. This indicated that, due to the influence of different factors such as imaging condition difference, land surface spatial heterogeneity, geometric registration and atmospheric correction, there were certain uncertainties in the reflectance data obtained by the UAV platform. The addition of reflectance noise improved the robustness of the two-layer matrix. In addition, the results also showed that the proposed method had better performance compared with PLSR and RF approaches.

4. Discussion

The remote sensing of LAI and Cab is challenging due to the cross-influence between them. The joint estimation approach based on a two-layer matrix proposed in this study can retrieve LAI and Cab simultaneously at the canopy level. Compared with the one-dimensional equations of individual VIs or VI ratios, matrix-based approaches can eliminate the influence of other parameters, such as leaf angle distribution, soil background, etc., as well as the cross-influence between LAI and Cab. The elimination of influencing is helpful for improving accuracy and reducing uncertainty in LAI and Cab. Compared with traditional physical modeling and machine learning methods, this approach is much simpler and more efficient.

Although the construction of vegetation indices has taken into account the influence of non-vegetation factors, such as soil background and atmosphere, it is still difficult to avoid the interaction between vegetation parameters. Therefore, in order to reduce the influence of different vegetation parameters during inversion, a variety of VI combinations are established. For example, a three-VI combination was established for calculating leaf mass per area (LMA), which takes into account the influence of leaf water content (Cw) in the SWIR region [41]. The TCARI–OSAVI combination was established for Cab retrieval from Sentinel-2 Data, which considers the influence of LAI on it [42]. However, the retrieval result of those VI combinations may still have uncertainty. In this paper, the retrieval result of LAI from UAV hyperspectral data using a single-layer matrix of OSAVI–MTCI had an R² of 0.55, which was still a low accuracy. The interaction between LAI and Cab was one of the factors affecting the inversion effect.

The estimation results of the two-layer matrix were obviously better than those of the single-layer matrix with and without noise. The reason why the two-layer matrix was superior to the single-layer matrix was that the two-layer matrix was combined with the most separable part of each single-layer matrix. This helped to reduce the cross-influence of LAI and Cab and improved the robustness of the inversion process. In addition, other vegetation and environmental parameters still affected the VI combinations and two-layer matrix. As shown in Section 2.2.3, parameters such as N, Cm and ALA had large effects on LAI and Cab estimation based on OSAVI–MTCI combination. A priori information can be used to limit the physical boundary of these parameters during simulation, so as to reduce the uncertainties [26].

Since the reflectance data are obtained after payload calibration, atmospheric correction and other processes, there is some uncertainty in the reflectance data. Therefore, the influence of reflectance measurement uncertainty on inversion was tested by adding or not adding noise tested in this study. If no noise was added during matrix modeling, the estimation accuracies of LAI and Cab retrieved from data with noise will be greatly reduced, as shown in Section 3.1. This showed that the uncertainty of the remotely sensed data can increase the uncertainty of retrieved LAI and Cab. If the phenomena of the same object with a different spectrum and the same spectrum with a different object were considered, noise was added in matrix modeling (namely the influence of LAI and Cab changes on reflectance). That is, the value in each cell of the matrix fluctuated within the noise range of 5% to generate mean value and standard deviation, which not only solved the ill-posed problem of physical model inversion, but also presented the uncertainty of the inversion.

The experiment for the joint estimation of the LAI and Cab from the simulated reflectance and UAV hyperspectral data both demonstrated that a two-layer matrix achieved better accuracies than the one-dimensional relationships of VIs and a single-layer matrix. Compared with machine learning methods such as RF and PLSR, the proposed method performed better, too. As a method similar to look-up table, the two-layer matrix can be used directly for LAI and Cab inversion without the need for prior measurements for training.

5. Conclusions

This paper proposed a novel approach for joint estimation of LAI and Cab using a two-layer matrix based on PROSAIL simulations. In this approach, LAI and Cab were directly obtained from a cell of the matrices, avoiding the possible distortion caused by the regression relationships of the VIs with LAI and Cab. This retrieval also reduced the cross-influence between LAI and Cab.

In this paper, the effectiveness of the approach for joint estimation of LAI and Cab was examined using simulated data and UAV hyperspectral data. Compared with the traditional single VIs and single-layer matrices, the two-layer matrix estimated LAI and Cab simultaneously, which improved inversion efficiency and accuracy. The results of adding noise showed that the simulation effect can be improved by considering the influence of noise when constructing the matrix. The accuracy of LAI and Cab retrieval from UAV hyperspectral data was significantly improved; the RMSE was reduced from 1.03 m²/m² to 0.91 m²/m² for LAI, and from 14.2 μg/cm² to 11.7 μg/cm² for Cab.