1.

Introduction

Grasslands cover about 40% of the earth’s terrestrial surface and thereby perform significant ecosystem functions, such as carbon sequestration,¹ and so are economically important.² A wide range of grasslands, such as intensive and extensive pastures and hay and silage meadows, are the basis of the world’s meat and milk production and of biofuel production and fibers.³ Of key importance in understanding the high spatio-temporal dynamics of cultivated and natural grasslands and adjusting management decisions is the monitoring of quantity and quality of above-ground biomass.⁴ Above-ground grass biomass is strongly correlated to canopy height.^5–8 Canopy height is, thus, an important parameter in management decisions related to grasslands such as grazing rotation or harvesting time.^8–10 Furthermore, precision agriculture applications, such as site-specific fertilizer applications, rely on plant parameter data with a high spatial and temporal resolution.^4,10 While traditional destructive field measurement techniques (e.g., clipping and weighing) are widely applied in daily farming practice, they are time-consuming, labor-intensive, and limited in characterizing the spatial variability of sward characteristics, an indicator of grassland biomass.^8,11 Since the 1960s, a range of commercially adapted applications have been developed to nondestructively measure grassland biomass. The simplest forms of handheld devices are pasture rulers (or sward sticks), which measure the uncompressed sward height. Rising plate or disk meters (RPMs) measure the compressed sward height by integrating sward height and density over a defined area using a weighted disk.^9,12 Several studies have proven that RPM-based compressed sward height is also a robust predictor for grassland biomass.^5,8,13

With technological development in recent decades, these devices have become more sophisticated. Devices such as the GrassHopper (TrueNorth Technologies, Shannon, Ireland) or GrassOmeter (Monford AG Systems Ltd., Dublin, Ireland) utilize ultrasonic distance sensors that can be mounted on a stick or boot. In addition to handheld devices, sensors such as the Pasture Meter (C-Dax Agricultural Solutions, Palmerston North, New Zealand) or the Pasture Reader (Naroaka Enterprises, Narracan, Australia) are mounted at the rear or front of a vehicle to measure plant height using a tunnel-like sensor equipped with an optical array or an ultrasonic measurement device, respectively. In addition to systems providing structural measures of the sward, biomass can also be assessed using optical sensing systems in the visible to near-infrared spectral region, such as the Yara N-Sensor (Yara ASA, Oslo, Norway) or the GreenSeeker (Trimble Inc., Sunnyvale, California, USA), both handheld or tractor-mounted. These sensors allow direct data logging and transfer and are linked with global navigation satellite systems to create yield maps.

Those sensing systems provide good estimates of sward height or biomass. However, some challenges still remain for precision grassland management or large-scale grassland ecosystem monitoring: (1) limited spatial coverage, especially for handheld devices, and therefore limitations in characterizing within-field spatial variability of the sward, (2) necessity for heavy technical equipment, (3) limited access to the field due to grazing animals or protected species, and, for the vehicle-mounted sensors, potential disturbances at higher frequency for repeated measurements, and finally (4) limitations of applicability depending on field conditions (e.g., slope and soil moisture).⁶

Remote sensing methods offer a potential for rapid and automated measurements of plant parameters, such as biomass, nitrogen or chlorophyll content, in high spatial and temporal resolution on a variety of spatial scales, especially for agricultural applications. These methods include digital imaging (hyperspectral, multispectral, RGB, radar), photogrammetry, laser scanning, and combinations of various sensors on different platforms.^14,15 Numerous studies have investigated using satellite remote sensing to estimate plant parameters. They derived crop biomass from spectral information of canopy reflectance.^14–19 However, most satellite systems with high spatial resolution ($<5\mathrm{m}$) are operated commercially, and thus the costs of image acquisition for short revisit times can become a constraining factor.²⁰

The rapid development of sensor and platform technology, especially in the field of unmanned aerial vehicles (UAVs), and of small high-resolution camera systems (standard RGB, multispectral, or hyperspectral) has opened up applications that support tasks in experimental fields and large farm areas, especially tasks related to precision grassland management and ecosystem monitoring.^9,21 In addition, the use of these technologies has been accelerated by a rapid development of user-friendly computer programs for three-dimensional scene reconstruction from aerial imagery by Structure-from-Motion (SfM) and Multi-View-Stereopsis (MVS) algorithms.²² Especially the flexible application of UAVs has benefitted agricultural monitoring approaches with high spatial and temporal resolution and multiple sensors.²³

In recent years, structural plant parameters, such as plant height, have become the focus of UAV-based remote sensing approaches for crop monitoring. Plant height derived from multitemporal canopy surface models (CSMs) has been studied as a robust estimator of biomass.^24–28 In addition, crop biomass has been estimated from spectral information from UAV-based standard RGB and multispectral or hyperspectral cameras.^29–34 Furthermore, some studies have investigated combining structural and spectral features from ground-based or UAV-based sensors to estimate crop parameters.^7,24,35–38

Those approaches have mainly been adopted to monitor arable crops, where spatial heterogeneity is often lower than in grasslands,¹⁰ which typically has a high spatio-temporal heterogeneity as a result of the considerably different floristic compositions and the co-occurrence of different phenological stages. This heterogeneity poses challenges in remote sensing-based estimation of biomass quantity and quality.³⁹ Some studies have faced these challenges by deploying ground-based measurements using ultrasonic sensors,^39,40 light detection and ranging,^41,42 or a combination of structural and spectral features.^{7,36,39,43,44} Furthermore, a few studies have been published on deploying structural features such as height or volume derived from UAV-based camera systems for grassland height or biomass estimation.^{6,28,36,45–48} However, to our knowledge, no study has utilized low-cost UAV-system-based sward height and vegetation indices (VIs) to estimate grassland biomass for at least one entire growing season with three consecutive growths in mid-latitude Europe.

Therefore, the main objectives of the present study are (1) to evaluate the potential of low-cost UAV-based CSMs to monitor sward height as an indicator of biomass, (2) to evaluate the potential of VIs from low-cost UAV-based RGB digital image data, and (3) to compare those methods with established methods for biomass monitoring such as RPMs and spectroradiometer-based narrowband VIs.

2.

Data and Methods

2.1.

Study Site

The study was conducted in 2017 on an experimental grassland site in Ersdorf (N 50°34′56.4″, E 6°59′21.1″) on the Campus Klein-Altendorf research facility of the Rheinische Friedrich Wilhelms University Bonn, Germany. The site is about 320 m above mean sea level and has a southeast exposition. The experimental site was established in March 2014 as a nitrogen fertilizer gradient experiment to develop a high range in biomass and plant height of perennial ryegrass (Lolium perenne). Further species are also present such as foxtail grass (Alopecurus sp.), timothy-grass (Phleum pratense) and clover (Trifolium sp.), and to a small amount dandelion (Taraxum officinale) and cocksfoot (Ranunculus sp.). The vegetation type was identified as a typical Lolio-Cynosuretum. The percentage contribution of species was limited so that dominance of species did not change significantly and the vegetation type remained the same.

The experimental setup comprised three growths with four cuts (including one equalization cut in early March) within one year, and the plots were set up in six fertilizer levels with three replicates, resulting in 162 plots of $1.5\times 3\mathrm{m}$ gross area (see Fig. 1). The plots were separated by 20-cm-wide border strips treated with herbicide. Calcium ammonium nitrate fertilizer (${\mathrm{CaH}}_4{\mathrm{N}}_4{\mathrm{O}}_9$, 27% N) was applied at the beginning of each growth. Thus, six fertilizer levels were established from 0 to $500\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-1}$ in increments of $100\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-1}$. Unfortunately, the, at-times, high activity of moles and the common vole disturbed parts of the sward significantly in several plots of replicate 1 (N1, N2, N5, N6) and replicate 2 (N4) (see Figs. 1 and 2).

Fig. 1

Location of the experimental site at the Campus Klein-Altendorf. The UAV-based orthomosaic in the detailed map of the experimental site was taken on the first sampling date of the first growth (April 26, 2017). Two exemplary heights of the GCPs show a slight slope. R, replicate; T1 to T3, sampling date 1 to 3 per growth and treatment; e.g. plot T1 in G2 was sampled on the first sampling date of the second growth.

Fig. 2

Close-up view (a) of swards with bent stalks due to vigorous growth in plots of replicate 2 treatment N5 (orthomosaic from August 22, 2017) and (b) of plots in replicate 1 with high rodent activity (orthomosaic from April 26, 2017).

The soil type was classified as a stagnosol (soil textural classes: silty clay loam/clay loam/clay). The annual precipitation in 2017 was 598.5 mm and the mean annual temperature was 10.8°C. Figure 3 displays the monthly mean precipitation and temperature for 2017 along with the long-term monthly means (1956 to 2015). April, May, and June of 2017 were hotter and dryer than the long-term averages.

Fig. 3

Monthly temperature (°C) and precipitation (mm) averages (a) for 2017 and (b) for the long-term average from 1956 to 2015. G1, G2, and G3 indicate duration of growths in 2017. Data are taken from the local weather station at Campus Klein-Altendorf.

Table 1 lists all sampling dates for the growing season in 2017. For each growth, three sampling dates were scheduled. Data collection included UAV campaigns, spectroradiometer measurements of canopy reflectance, and RPM measurements of compressed sward height and destructive biomass sampling.

Table 1

Sampling dates for UAV campaigns, reflectance measurements, and reference field data.

Growth	Sampling date number	Date	UAV	ASD	RPM	Biomass
G1	T0	2017-03-29	x	—	—	—
	T1	2017-04-26	x	—	x	x
	T2	2017-05-03	x	x R. 2 and 3	x	x
	T3	2017-05-10	x	x R. 2 and 3	x	x
G2	T1	2017-06-07	x	x	x	x
	T2	2017-06-20	x	x	x	x
	T3	2017-06-28	x	—	x	x
G3	T0	2017-07-04	x	—	—	—
	T1	2017-08-02	x	x	x	x
	T2	2017-08-09	x	x R. 2 and 3	x	x
	T3	2017-08-22	x	x	x	x

RPM, rising plate meter; ASD, ASD FieldSpec3; R: replicate; x/—: data/no data; For G1T2, G1T3, and G3T2 only replicates 2 and 3 were measured with the ASD instrument.

3.

Results

3.1.

Orthomosaics and Canopy Surface Models

The orthomosaics and CSMs for the three sampling dates of the first growth are presented in Fig. 4. Both display a good visual correlation to the fertilizer treatments and a detailed representation of sward height is achieved with high spatial resolution. The disturbances by rodents as mentioned above became visible in the CSMs and orthomosaics. In the orthomosaics of Fig. 4, the effect of a rainfall event on G1-T2 shortly before image acquisition can be seen: the wet soil appears darker than the dry soil on the first and third sampling dates.

Fig. 4

CSMs (left panels) and orthomosaics (right panels) for the first growth in 2017: (a) and (b) T1 April 26, 2017; (c) and (d) T2 May 3, 2017; (e) and (f) T3 May 10, 2017. Lower sward heights from previously sampled plots or from plots affected by high rodent activity are clearly distinguishable.

3.2.

Reference Measurements

Due to vigorous growth in G3, grass stalks were slightly lodging in the highest fertilizer treatments (N5 and N6) on G3-T3 (see Fig. 2). Plots with high rodent activity were excluded from the analysis, since they produced outliers with unreasonable biomass values.

Table 5 summarizes the descriptive statistics of the biomass measurements. G1 and G3 yielded the highest biomass values, while G2 was affected by drought and thus had lower tiller density and biomass development.

Table 5

Descriptive statistics of the biomass measurements for each sampling date per growth. G, growth; T, sampling date number; min, minimum; max, maximum; sd, standard deviation.

		Date	DBM t ha−1				FBM t ha−1
			min	max	mean	sd	min	max	mean	sd
G1	T1	2017-04-26	1.33	2.87	2.02	0.52	4.49	10.64	7.10	2.12
	T2	2017-05-03	1.45	3.83	2.44	0.69	5.01	15.94	9.41	3.22
	T3	2017-05-10	1.41	4.37	2.74	0.97	4.56	15.60	9.35	3.60
G2	T1	2017-06-07	0.58	2.27	1.37	0.47	1.95	7.56	4.46	1.66
	T2	2017-20-06	0.55	1.70	1.00	0.34	1.90	5.36	3.22	1.03
	T3	2017-06-28	1.16	2.36	1.74	0.40	4.29	9.83	6.65	1.66
G3	T1	2017-08-02	0.68	2.17	1.43	0.39	2.21	9.85	5.81	2.24
	T2	2017-08-09	1.67	3.57	2.34	0.46	5.30	13.22	8.03	2.00
	T3	2017-08-22	2.15	4.10	3.32	0.60	6.94	16.40	11.91	2.68

Table 6 summarizes the descriptive statistics for the RPM sward height reference measurements (${\mathrm{SH}}_{\mathrm{RPM}}$) and the mean sward height derived from CSM (${\mathrm{SH}}_{\mathrm{CSM}}$). The ${\mathrm{SH}}_{\mathrm{CSM}}$ resulted in a higher range of values and a higher standard deviation than the ${\mathrm{SH}}_{\mathrm{RPM}}$.

Table 6

Descriptive statistics of the sward height reference measurements (RPM) and mean sward height from CSMs for each sampling date per growth. G, growth; T, sampling date number; min, minimum; max, maximum; sd, standard deviation.

		Date	Sward height RPM (cm)				Mean sward height CSM (cm)
			min	max	mean	sd	min	max	mean	sd
G1	T1	2017-04-26	8.14	13.42	10.72	1.67	2.47	14.00	9.43	3.17
	T2	2017-05-03	8.90	18.26	12.29	2.59	10.97	22.52	16.14	3.22
	T3	2017-05-10	9.84	19.10	14.87	3.34	6.20	22.41	12.88	5.12
G2	T1	2017-06-07	5.82	11.20	8.42	1.46	8.43	19.10	15.39	3.50
	T2	2017-20-06	6.18	9.36	7.27	0.92	4.74	13.66	8.31	2.84
	T3	2017-06-28	6.44	9.74	7.91	0.96	14.94	22.80	18.79	2.68
G3	T1	2017-08-02	7.40	12.96	9.41	1.52	2.38	13.25	7.27	2.70
	T2	2017-08-09	8.18	13.84	10.73	1.83	10.81	18.98	14.26	2.27
	T3	2017-08-22	9.74	16.28	13.92	1.89	9.43	21.58	14.52	3.08

Table 7 lists the PCCs for the CSM-based height features ${\mathrm{SH}}_{\mathrm{p}90}$ and ${\mathrm{SH}}_{\text{mean}}$ and the RPM reference measurements of DBM and FBM. Since the ${\mathrm{SH}}_{\mathrm{p}90}$ correlated slightly better than ${\mathrm{SH}}_{\text{mean}}$ for DBM and FBM, only the ${\mathrm{SH}}_{\mathrm{p}90}$ is described in the following section. For completeness, a sensitivity analysis of all single features by growths, by growths and treatment, and by growths and sampling date is included in the Appendix (Tables 12–19).

Table 7

PCCs for the CSM-derived height features (SHmean: CSM-based mean sward height, SHp90: CSM-based 90th percentile sward height) and rising plate meter (SHRPM) reference measurements to DBM and FBM for all growths.

	Growth 1		Growth 2		Growth 3
	DBM	FBM	DBM	FBM	DBM	FBM
${\mathrm{SH}}_{\mathrm{RPM}}$	0.94	0.89	0.78	0.68	0.89	0.92
${\mathrm{SH}}_{\text{mean}}$	0.80	0.85	0.81	0.79	0.75	0.65
${\mathrm{SH}}_{\mathrm{p}90}$	0.87	0.90	0.83	0.83	0.78	0.68

P<0.001 for all correlations.

The relationships between DBM and ${\mathrm{SH}}_{\mathrm{RPM}}$ and ${\mathrm{SH}}_{\mathrm{p}90}$ are displayed in Fig. 5. ${\mathrm{SH}}_{\mathrm{RPM}}$ shows a strong linear relationship to DBM with little deviation from the regression line for all three growths.

Fig. 5

Scatterplots of (a) mean CSM-based sward height versus DBM yield, (b) 90th percentile CSM-based sward height versus DBM yield, (c) mean RPM-based sward height versus DBM yield, and (d) 90th percentile CSM-based sward height versus mean RPM-based sward height. $P<0.001$ for all correlations.

The ${\mathrm{SH}}_{\mathrm{p}90}$ feature shows larger deviation from the regression line, while the second growth falls outside the pattern of the first and third growths. Nevertheless, the $R^2$ for the second growth was higher for ${\mathrm{SH}}_{\mathrm{p}90}$, compared to the ${\mathrm{SH}}_{\mathrm{RPM}}$. ${\mathrm{SH}}_{\mathrm{RPM}}$ and ${\mathrm{SH}}_{\mathrm{p}90}$ were only moderately well correlated.

Figure 6 presents the relationship between selected VIs and DBM. The UAV-based NGRDI performed moderately well for the first growth with an $R^2$ of 0.58 but shows no correlation for the second and third growths. In comparison, the ASD-based nNGRDI performed more consistently for the three growths (see Table 15). Slight lodging of the grass sward in six plots with high fertilizer treatments at G3-T3 led to a higher deviation from the regression line. The spectroradiometer-based NDVI performed moderately well. In comparison to the NGRDI, the ExGI correlated better with DBM in G3 since the slightly lodging sward in six plots of treatments N5 and N6 seemed to have only a minor effect on the ExGI values.

Fig. 6

Scatterplots of selected VIs for all growths: (a) NDVI versus DBM yield ($P<0.001$ for all growths), (b) ExGI versus DBM yield ($P<0.001$, 0.208, and $<0.001$ for G1, G2, and G3, respectively), (c) ASD-based narrowband NGRDI versus DBM yield ($P<0.001$ for all growths), and (d) orthomosaic-based NGRDI versus DBM yield ($P<0.001$, 0.067, 0.521 for G1, G2, and G3, respectively).

3.3.

Cross-Validation Results of Simple Linear Regression for Estimation of Dry Biomass and Fresh Biomass

The predictive accuracy of the height features and VIs for biomass estimation was assessed using LOO-CV. Table 8 displays the cross-validation results of the bivariate regression models of DBM and each feature by growth. The ${\mathrm{SH}}_{\mathrm{RPM}}$ performed best for G1 and G3 in estimating DBM, with an $R^2$ of 0.87 (RMSE: $0.274\mathrm{t}{\mathrm{ha}}^{-1}$) and 0.78 (RMSE: $0.416\mathrm{t}{\mathrm{ha}}^{-1}$), respectively, while the ${\mathrm{SH}}_{\mathrm{p}90}$ performed best for G2 with an $R^2$ of 0.66 (RMSE: $0.287\mathrm{t}{\mathrm{ha}}^{-1}$). The performance of the UAV-based ${\mathrm{VI}}_{\mathrm{RGB}}$ depended on the growth. NGRDI and VARI showed a moderate performance with an $R^2$ of 0.54 and 0.58, respectively, for G1 but showed no correlation for G2 and G3. Excluding G3-T3 (slightly lodging sward in six plots of treatments N5 and N6) led to a better performance of both indices with an $R^2$ of 0.42 (RMSE: $0.469\mathrm{t}{\mathrm{ha}}^{-1}$) and 0.56 (RMSE: $0.409\mathrm{t}{\mathrm{ha}}^{-1}$), respectively. Similarly, the $R^2$ of the RGBVI increased to 0.55 (RMSE: $0.416\mathrm{t}{\mathrm{ha}}^{-1}$). However, the RGBVI showed no predictive accuracy for G1 and G2. The ExGI had little-to-no predictive accuracy for G1 and G2 but an $R^2$ of 0.55 (RMSE: $0.601\mathrm{t}{\mathrm{ha}}^{-1}$) for G3. The ASD-based narrowband ${\mathrm{VI}}_{\mathrm{RGB}}$ performed well for DBM estimation for G1 and G3, except for the ExGI. However, only a weak correlation was observed for G2.

Table 8

Statistics of cross-validation results for bivariate linear regression of DBM yield against each predictor variable per growth. R2CV, Cross-validation coefficient of determination; RMSECV, cross-validation root mean squared error; n, number of samples.

Estimator		Growth 1		Growth 2		Growth 3
		R2CV	RMSECV (t ha−1)	R2CV	RMSECV (t ha−1)	R2CV	RMSECV (t ha−1)
	RPM	0.87	0.274	0.58	0.320	0.78	0.416
UAV	${\mathrm{SH}}_{\text{mean}}$	0.61	0.475	0.62	0.305	0.53	0.612
	${\mathrm{SH}}_{\mathrm{p}90}$	0.73	0.394	0.66	0.287	0.57	0.586
	RGBVI	0.06	0.740	0.06	0.513	0.32	0.739
	NGRDI	0.54	0.518	0.02	0.490	0.11	0.927
	VARI	0.58	0.494	0.00	0.502	0.09	0.924
	ExGI	0.24	0.664	0.01	0.508	0.55	0.601
	GrassI	0.61	0.475	0.61	0.306	0.54	0.610
	${\mathrm{ExGI}+\mathrm{SH}}_{\mathrm{p}90}$	0.73	0.394	0.66	0.287	0.58	0.582
	$n$	43		42		45
ASD	NDVI	0.83	0.360	0.55	0.295	0.32	0.767
	OSAVI	0.68	0.490	0.51	0.306	0.39	0.723
	RDVI	0.56	0.577	0.47	0.321	0.40	0.719
	REIP	0.86	0.322	0.29	0.369	0.02	0.924
	NDREI	0.87	0.316	0.56	0.291	0.04	0.915
	nRGBVI	0.64	0.520	0.28	0.373	0.58	0.602
	nNGRDI	0.76	0.422	0.37	0.349	0.57	0.605
	nVARI	0.78	0.404	0.38	0.345	0.49	0.664
	nExGI	0.44	0.648	0.01	0.459	0.03	0.923
	$n$	22		30		41

The composite indices GrassI and $\mathrm{ExGI}+{\mathrm{SH}}_{\mathrm{p}90}$ performed as expected in the range of the best-performing component (height features from CSM). The ${\mathrm{VI}}_{\mathrm{VNIR}}$ performed best for G1. The NDVI, REIP, and NDREI had an $R^2$ of 0.83, 0.86, and 0.87, respectively, while the RDVI performed weak for G1. For G2, the ${\mathrm{VI}}_{\mathrm{VNIR}}$ showed moderate predictive accuracy, while the REIP only had an $R^2$ of 0.29. In G3, the red-edge-based indices REIP and NDREI showed no correlation, and the NDVI, OSAVI, and RDVI showed only a weak predictive accuracy. Similar to the ${\mathrm{VI}}_{\mathrm{RGB}}$, the exclusion of G3-T3 increased the performance of the NDVI, REIP, and NDREI, although only a low predictive accuracy could be achieved (NDVI: $R^2$ 0.39, RMSE $0.512\mathrm{t}{\mathrm{ha}}^{-1}$; REIP: $R^2$ 0.24, RMSE $0.574\mathrm{t}{\mathrm{ha}}^{-1}$; NDREI $R^2$ 0.26, RMSE $0.567\mathrm{t}{\mathrm{ha}}^{-1}$).

Table 9 presents the cross-validation results of the FBM prediction based on simple linear regression. As for the DBM prediction, ${\mathrm{SH}}_{\mathrm{RPM}}$ is the best predictor for FBM for G1 and G3 with an $R^2$ of 0.78 (RMSE: $1.448\mathrm{t}{\mathrm{ha}}^{-1}$) and 0.83 (RMSE: $1.367\mathrm{t}{\mathrm{ha}}^{-1}$), respectively. However, ${\mathrm{SH}}_{\mathrm{p}90}$ has a slightly higher $R^2$ of 0.79 (RMSE: $1.421\mathrm{t}{\mathrm{ha}}^{-1}$) for G1 and outperforms ${\mathrm{SH}}_{\mathrm{RPM}}$ for G2 with an $R^2$ of 0.65 (RMSE: $1.166\mathrm{t}{\mathrm{ha}}^{-1}$).

Table 9

Statistics of cross-validation results for bivariate linear regression of FBM against each predictor variable per growth. R2CV, cross-validation coefficient of determination; RMSECV, cross-validation root mean squared error; n, number of samples.

Estimator		Growth 1		Growth 2		Growth 3
		R2CV	RMSECV (t ha−1)	R2CV	RMSECV (t ha−1)	R2CV	RMSECV (t ha−1)
	RPM	0.78	1.448	0.42	1.504	0.83	1.367
UAV	${\mathrm{SH}}_{\text{mean}}$	0.71	1.670	0.60	1.253	0.38	2.633
	${\mathrm{SH}}_{\mathrm{p}90}$	0.79	1.421	0.65	1.166	0.43	2.526
	RGBVI	0.12	2.897	0.55	2.070	0.31	2.773
	NGRDI	0.68	1.756	0.06	1.926	0.00	3.397
	VARI	0.71	1.648	0.02	1.977	0.01	3.418
	ExGI	0.35	2.491	0.12	2.060	0.51	2.331
	GrassI	0.71	1.667	0.59	1.261	0.38	2.624
	${\mathrm{ExGISH}}_{\mathrm{p}90}$	0.79	1.410	0.65	1.169	0.43	2.513
	$n$	43		42		45
ASD	NDVI	0.84	1.393	0.65	0.868	0.43	2.576
	OSAVI	0.81	1.515	0.67	0.842	0.57	2.252
	RDVI	0.73	1.812	0.65	0.868	0.59	2.203
	REIP	0.89	1.162	0.36	1.180	0.12	3.213
	NDREI	0.89	1.172	0.67	0.842	0.15	3.161
	nRGBVI	0.59	2.239	0.36	1.188	0.53	2.354
	nNGRDI	0.73	1.801	0.47	1.074	0.61	2.136
	nVARI	0.77	1.656	0.49	1.053	0.57	2.237
	nExGI	0.34	2.839	0.00	1.525	0.00	3.474
	$n$	22		30		41

The best-performing ${\mathrm{VI}}_{\mathrm{RGB}}$ for G1 were the NGRDI and VARI, while both indices obviously failed to estimate FBM for G2 and G3. The RGBVI showed a better predictive accuracy for FBM ($R^2$: 0.55, RMSE: $2.070\mathrm{t}{\mathrm{ha}}^{-1}$) than for DBM in G2, while the ExGI followed a similar pattern for FBM prediction as for DBM, with the best performance in G3 ($R^2$: 0.51, RMSE: $2.331\mathrm{t}{\mathrm{ha}}^{\mathrm{v}1}$). The performance of the ASD-based narrowband ${\mathrm{VI}}_{\mathrm{RGB}}$ for FBM estimation was similar to that for DBM estimation. As expected, the composite indices GrassI and ExGI+SHp90 performed in the range of the best-performing single feature (height features from CSM). The ${\mathrm{VI}}_{\mathrm{VNIR}}$ performed best for G1. All indices had a high $R^2$ value above 0.70, and both the red-edge-based indices REIP and NDREI had an $R^2$ of 0.89. Also, for G2, the $R^2$ values for FBM prediction of the ${\mathrm{VI}}_{\mathrm{VNIR}}$ were 0.65 or slightly higher, except for the REIP, which only had an $R^2$ of 0.36 (RMSE: $1.180\mathrm{t}{\mathrm{ha}}^{-1}$). For G3, the best-performing indices were the OSAVI and RDVI with an $R^2$ of 0.57 (RMSE: $2.252\mathrm{t}{\mathrm{ha}}^{-1}$) and 0.59 (RMSE: $2.203\mathrm{t}{\mathrm{ha}}^{-1}$), respectively. The lowest RMSE values were observed for the ${\mathrm{VI}}_{\mathrm{VNIR}}$ for G2, where all indices except the REIP had an RMSE below $0.900\mathrm{t}{\mathrm{ha}}^{-1}$.

Observed and predicted DBM for selected features are displayed in Fig. 7. ${\mathrm{SH}}_{\mathrm{RPM}}$ is closest to the 1:1 line, while predictions based on NDVI and ${\mathrm{SH}}_{\mathrm{p}90}$ deviate more from the regression line, especially for higher DBM values.

Fig. 7

Scatterplots of observed versus predicted DBM yield for (a) ${\mathrm{SH}}_{\mathrm{p}90}$ as predictor variable, (b) NDVI as predictor variable, and (c) ${\mathrm{SH}}_{\mathrm{RPM}}$ as predictor variable. Cross-validation $R^2$ values are presented.

4.

Discussion

The primary aims of this study were (1) evaluating the ability of height features and VI features derived from UAV-based image data to predict DBM and FBM in grassland and (2) comparing the performance of these features with established VIs from the VNIR spectral region and RPM reference measurements. Although Bareth and Schellberg⁶ found that ${\mathrm{SH}}_{\mathrm{RPM}}$ and the ${\mathrm{SH}}_{\mathrm{CSM}}$ correlated well, in this study, both features only correlate at a low-to-medium level depending on the growths. The ${\mathrm{SH}}_{\mathrm{CSM}}$ features tend to overestimate the manually measured ${\mathrm{SH}}_{\mathrm{RPM}}$, which was expected due to the compression of the sward when using the RPM.⁶ Furthermore, the choice of the RPM instrument might have influenced the correlation due to a different disk weight and diameter.

Recent studies using SfM-MVS to derive canopy height models for grassland have obtained reference measurements in the field with a height stick or a ruler. Dependent on the sward structure, species composition, and growing stage, Grüner et al.⁴⁶ achieved $R^2$ values of 0.56 and 0.70 and Viljanen et al.⁷ report $R^2$ values of 0.61 to 0.93 when comparing ${\mathrm{SH}}_{\mathrm{CSM}}$ with manual reference measurements from height sticks. Forsmoo et al.⁴⁷ achieved $R^2$ values of 0.64 and 0.72 when comparing reference measurements from RPMs and ${\mathrm{SH}}_{\mathrm{CSM}}$ for a grassland field of mainly perennial ryegrass and clover in southwest England. Higher correlation between RPM reference measurements and low-cost ${\mathrm{SH}}_{\mathrm{CSM}}$ ($R^2$ of 0.83 to 0.91) were reported by Bareth and Schellberg⁶ for three consecutive years in the long-term Rengen Grassland Experiment (RGE) in Germany. In contrast to manual measurements, the CSM approach captures the spatial variability of the plant height in much finer detail.

To obtain sward height, the CSM-approach described by Bendig et al.⁴⁹ was used. Table 2 shows that the errors of the reconstructed DSMs in $X$, $Y$, and $Z$ directions were within the range of up to 7.3 cm, which clearly have an impact on calculating sward height when $\mathrm{DSM}\text{-}{\mathrm{T}}_0$ is subtracted from $\mathrm{DSM}\text{-}{\mathrm{T}}_{1\text{-}n}$. The high error rates may be due to several factors such as the suboptimal placement of GCPs or the sensor geometry of the DJI P4A camera. Owing to overgrowing of the border strips between the plots, it was not possible to classify a DTM representing the ground for every sampling date. For small-scale field experiments, DTM classification or measuring ground surface points with an real time kinematic differential GPS may be more feasible. Nevertheless, the CSM approach by Bendig et al.⁴⁹ may be more applicable to actual field conditions since ground surface classification may not be viable for a fully covered grassland field of several hectares with varying topography. However, the CSM approach by Bendig et al⁴⁹ requires that the data for the base model are acquired over the bare soil surface of the field, which is not feasible for grassland fields due to permanent vegetation coverage.

The high rodent activity in a few plots per growth led to higher uncertainties in ${\mathrm{SH}}_{\mathrm{CSM}}$ calculation, which, in some cases, resulted in negative sward height values and unreasonable biomass values. The decision to exclude these plots improved the biomass estimation, although these disturbances are likely to occur in the field.

Biomass prediction models based on LOO-CV results showed moderate accuracies depending on the growths when using the height features (${\mathrm{SH}}_{\text{mean}}$ and ${\mathrm{SH}}_{\mathrm{p}90}$). Viljanen et al.⁷ reported correlation coefficients between 0.75 and 0.98 for biomass prediction using height features from CSM but concluded that biomass prediction accuracy was dependent on the density and growth stage of the sward. Näsi et al.³⁶ reported correlation coefficients for grassland biomass prediction between 0.10 and 0.41. Grüner et al.⁴⁶ achieved $R^2$ values of 0.46 and up to 0.87 depending on the sward composition for mixed legume-grass swards and pure legumes and grass stands. Roth and Streit⁵⁰ analyzed different cover crops, and the prediction accuracy for biomass improved when lodging plants were excluded from the analysis. Our findings are, thus, comparable to other studies. Since G2 was affected by drought (see Fig. 3), low biomass, low tiller density, and low sward heights in G2 led to a weak correlation when all growths were combined in a single regression analysis (see Fig. 8).

Fig. 8

Relationship between (a) 90th percentile CSM-based sward height and DBM yield for all growths in 2017 and (b) only for growths 1 and 2.

Similar results were also observed by Grüner et al.⁴⁶ when pure grass stands from the second growth in 2017 were excluded. These findings indicate that challenges remain when using SfM-MVS reconstruction to determine a global model for biomass estimation in pure grass swards.

The varying accuracy especially of the ${\mathrm{VI}}_{\mathrm{RGB}}$ from the orthomosaics as biomass predictors might be explained by several factors such as the varying soil color due to rainfall events or drought, bent stalks or flowers. Owing to vigorous growth in plots of treatments N5 and N6 in G3-T3, stalks tended to bend over, and as a result, a higher proportion of the spectral signature of the stalks was captured instead of the desired nadir view of the canopy (see Fig. 2). In this study, the correlation of the UAV-based ${\mathrm{VI}}_{\mathrm{RGB}}$ by growth and treatment did not yield more robust results than the analysis of the features excluding the treatment effects. Furthermore, no clear pattern between treatments and VIs is recognizable in the sensitivity analysis (see Tables 14–16). For clarification, an extended sensitivity analysis of all single features can be found in the Appendix.

Dependent on the sampling date, the illumination conditions were stable but either sunny or overcast, which probably contributed significantly to the UAV-based ${\mathrm{VI}}_{\mathrm{RGB}}$ values. This observation was also shown by Rasmussen et al.⁷¹ for wheat. Furthermore, the effect of rainfall shortly before image acquisition probably influenced the ${\mathrm{VI}}_{\mathrm{RGB}}$ values due to changing soil color (see Fig. 1) The study by Viljanen et al.⁷ achieved correlation coefficients of 0.70 to 0.85 for the correlation of biomass and RGBVI, NGRDI, and ExGI. However, for the random forest classifier, the ${\mathrm{VI}}_{\mathrm{RGB}}$ was not among the most important features in estimating grass biomass.

The choice of multiple linear regression (one VI feature + one height feature) was motivated by the search for a robust and simple model that can be applied to different growths per season. Initial tests on different interaction terms for the MLR (such as VI feature × height feature) did not yield more robust results.

The performance of the composite VIs GrassI and $\mathrm{ExGI}+{\mathrm{SH}}_{\mathrm{p}90}$ did not yield significantly better results than the multivariate regression of the respective ${\mathrm{VI}}_{\mathrm{RGB}}$ and height features, as they rely on the quality of both features. Owing to the different viewing geometry and technique of the ASD-based VIs, these features are not directly comparable to the results of the UAV-based ${\mathrm{VI}}_{\mathrm{RGB}}$. As expected, VIs from a well-calibrated instrument, and including the NIR spectral region, performed better and more consistently over all growths than the UAV-based ${\mathrm{VI}}_{\mathrm{RGB}}$.^24,38 However, all ASD-based VIs had a decreasing accuracy from G1 to G3 and were affected by slightly lodging sward in six plots with higher N-treatments, as mentioned already. The ASD-based narrowband ${\mathrm{VI}}_{\mathrm{RGB}}$ validates the respective VIs as viable for biomass prediction in grassland and demonstrate the challenges associated with uncalibrated sensors such as the camera used in this study.

The heterogeneous sward structure with high spatiotemporal variability compared to crops led to varying performance for biomass estimation depending on the growths and choice of predictor variable. It is therefore essential to study this topic further and to evaluate different swards under varying conditions and sites and over multiple years. A promising approach is the application of high-resolution multispectral imaging sensors, such as the MicaSense RedEdge camera (Micasense Inc., Seattle), which can readily be employed on a UAV. Furthermore, comparing simple regression techniques with more sophisticated algorithms such as random forest should be explored for larger datasets of grassland biomass and UAV-based structural and spectral features.

5.

Conclusions

Estimating biomass in high spatial and temporal resolution is a key component in precision agriculture applications and ecosystem monitoring. In this study, SfM-MVS-derived features of sward height and VIs from UAV-based images were assessed as predictors of grassland above-ground biomass and compared to established narrowband VIs from spectroradiometer measurements. The application of UAV-based imaging sensors serves as a fast and nondestructive method for data acquisition with high spatial and temporal resolution.

This study has shown that, especially for the first growth, which is considered the most important from an agronomical point of view, grassland biomass estimation by SfM-MVS-derived sward height from UAV-based images is feasible and provides an alternative means to manual measurement techniques such as RPMs or clipping. However, the results are influenced by various factors such as growing stage, sward composition, and biotic and abiotic factors, which need to be further investigated. Further research should be focused on integrating structural and spectral features for grassland biomass estimation and on generalizing models to different sites and years.

6.

Appendix

The appendix contains Tables 12–19 for sensitivity analysis (Pearson’s correlation coefficient, PCC) of the features (vegetation indices and canopy height) tested in this study for estimation of dry biomass yield (DBM).

Table 12

PCCs for DBM and UAV-based VIs by growth (G). ExGI, ExcessGreen Index; NGRDI, Normalized Green Red Difference Index; RGBVI, Red-Green-Blue Vegetation Index; VARI, Visible Atmospherically Resistant Index.

Table 13

PCCs for DBM and narrowband ASD-based VIs by growth (G). NDREI, Normalized Difference Red-Edge Index; NDVI, Normalized Difference Vegetation Index; OSAVI, Optimized Soil Adjusted Vegetation Index; RDVI, Renormalized Difference Vegetation Index; REIP, Red Edge Inflection Point; nExGI, ExcessGreen Index; nNGRDI, Normalized Green Red Difference Index; nRGBVI, Red-Green-Blue Vegetation Index; and nVARI, Visible Atmospherically Resistant Index.

Table 14

PCCs for DBM and UAV-based VIs and height features by growth (G) and treatment (N: kg N / ha). ExGI, ExcessGreen Index; NGRDI, Normalized Green Red Difference Index; RGBVI, Red-Green-Blue Vegetation Index; VARI, Visible Atmospherically Resistant Index; SHmean, mean sward height from CSMs; SHp90, 90th percentile sward height from CSM.

Table 15

PCCs for DBM and narrowband ASD-based RGB VIs by growth (G) and treatment (N: kg N / ha). nExGI, ExcessGreen Index; nNGRDI, Normalized Green Red Difference Index; nRGBVI, Red-Green-Blue Vegetation Index; nVARI, Visible Atmospherically Resistant Index.

Table 16

PCCs for DBM and ASD-based VNIR VIs by growth (G) and treatment (N: kg N / ha). NDREI, Normalized Difference Red-Edge Index; NDVI, Normalized Difference Vegetation Index; OSAVI, Optimized Soil Adjusted Vegetation Index; RDVI, Renormalized Difference Vegetation Index; REIP, Red Edge Inflection Point.

Table 17

PCCs for DBM and UAV-based VIs and height features by growth (G) and sampling date (T). ExGI, ExcessGreen Index; NGRDI, Normalized Green Red Difference Index; RGBVI, Red-Green-Blue Vegetation Index; VARI, Visible Atmospherically Resistant Index; SHmean, mean sward height from CSMs; SHp90, 90th percentile sward height from CSM.

Table 18

PCCs for DBM and narrowband ASD-based RGB VIs by growth (G) and sampling date (T). nExGI, ExcessGreen Index; nNGRDI, Normalized Green Red Difference Index; nRGBVI, Red-Green-Blue Vegetation Index; nVARI, Visible Atmospherically Resistant Index.

Table 19

PCCs for DBM and narrowband ASD-based VNIR VIs by growth (G) and sampling date (T). NDREI, Normalized Difference Red-Edge Index; NDVI, Normalized Difference Vegetation Index; OSAVI, Optimized Soil Adjusted Vegetation Index; RDVI, Renormalized Difference Vegetation Index; REIP, Red Edge Inflection Point.