Fusion of Plant Height and Vegetation Indices for the Estimation of Barley Biomass

Abstract

Plant biomass is an important parameter for crop management and yield estimation. However, since biomass cannot be determined non-destructively, other plant parameters are used for estimations. In this study, plant height and hyperspectral data were used for barley biomass estimations with bivariate and multivariate models. During three consecutive growing seasons a terrestrial laser scanner was used to establish crop surface models for a pixel-wise calculation of plant height and manual measurements of plant height confirmed the results (R² up to 0.98). Hyperspectral reflectance measurements were conducted with a field spectrometer and used for calculating six vegetation indices (VIs), which have been found to be related to biomass and LAI: GnyLi, NDVI, NRI, RDVI, REIP, and RGBVI. Furthermore, biomass samples were destructively taken on almost the same dates. Linear and exponential biomass regression models (BRMs) were established for evaluating plant height and VIs as estimators of fresh and dry biomass. Each BRM was established for the whole observed period and pre-anthesis, which is important for management decisions. Bivariate BRMs supported plant height as a strong estimator (R² up to 0.85), whereas BRMs based on individual VIs showed varying performances (R²: 0.07–0.87). Fused approaches, where plant height and one VI were used for establishing multivariate BRMs, yielded improvements in some cases (R² up to 0.89). Overall, this study reveals the potential of remotely-sensed plant parameters for estimations of barley biomass. Moreover, it is a first step towards the fusion of 3D spatial and spectral measurements for improving non-destructive biomass estimations.

1. Introduction

Over the past several decades remote sensing has increased in importance for precision agriculture [1,2,3]. Since the world population is expected to increase by more than one third until 2050 a main goal is shrinking the gap between potential and current yield [4,5]. Field management strategies in precision agriculture that aim to maximize yield must involve a reasonable use of natural resources and have to take spatial and temporal variabilities into account [6], as agricultural production is influenced by the physical landscape, climatic variables, and agricultural management practices [2]. Studies reveal that grain yield is correlated with total biomass [7,8]. A quantitative measure is the harvest index, which expresses yield vs. total biomass [9]. Moreover, adequate crop condition in early growing stages could buffer the yield against environmental stresses, such as droughts, during later stages [10]. In-season, the nitrogen nutrition index, the ratio between actual and critical nitrogen (N) content, is widely used as a measure of the plant status [11]. The critical value is defined by a crop-specific N dilution curve, showing the relation between N concentration and biomass. Hence, an exact in-season acquisition of biomass is important in precision agriculture.

Since plant biomass cannot be determined non-destructively, other plant parameters are used as estimators. Therefore, remote sensing measurements enable an objective and accurate acquisition in a high temporal frequency [2]. A review of remote sensing methods for assessing biomass is given by Ahamed et al. [12]. At the field level, ground-based methods are commonly used to achieve sufficiently high resolutions and over the last several decades, several studies investigated the relationship between spectral reflectance measurements and crop characteristics. For extracting information, various vegetation indices (VIs) were developed from the reflectance in determined wavelengths. Two band VIs like the normalized difference vegetation index (NDVI) were traditionally used with multispectral broad band systems to estimate biomass or biomass related parameters, like LAI. Such VIs have been adapted to narrow band hyperspectral data and other band combinations [13,14,15,16]. Additionally, other VIs with more than two bands, such as the GnyLi, have been developed for the same purpose [17].

Moreover, active sensors based on light detection and ranging (LiDAR) have been increasingly used in vegetation studies since the 1980s [18]. Indeed, a main benefit of LiDAR is the very high resolution, which enables the acquisition of complex canopies [19]. In agricultural applications, for example, ground-based LiDAR methods, also known as terrestrial laser scanning (TLS), reveal potential for assessing plant height [20], leaf area index [21], crop density [22,23,24], or post-harvest growth [25]. Furthermore, the potential for estimating biomass with TLS is supported through studies on small grain cereals [26,27,28], sagebrush [29], and paddy rice [30,31]. The 3D architecture of single plants was modeled under laboratory conditions [32,33], however the transferability of those laboratory results to field conditions has not yet been shown.

Generally, the accuracy of estimations is a major issue, with the accuracy being limited when calculations are based on one estimator. Whilst biomass estimations based on VIs are affected by saturation effects [13,34,35], plant height may reach limitations when differences in plant height are low. Consequently, the fusion of multiple parameters should be examined to enhance estimations. So far, studies on the fusion of spectral and non-spectral information have been applied for characterization of forest ecosystems [36] and modeling of corn yield [37]. As both studies applied airborne methods, the spatial resolution was low. A ground-based multi-sensor approach for predicting biomass of grassland based on measurements of plant height, leaf area index (LAI), and spectral reflectance showed that combining multiple sensors can improve the estimation [34]. However, in that study, spectral data were not well-suited. Recently, the potential of the combined use of spectral and non-spectral ground-based measurements for estimating biomass was demonstrated for rice, maize, cotton, and alfalfa [15].

The overall aim of this study was to compare the potential of plant height (PH), VIs, and a fusion of PH and VIs for estimations of above ground fresh and dry barley biomass. More specifically, this study compares the potential of 3D spatial and spectral information for different time frames during the growing season and investigates if a fusion of both can improve the estimation. Therefore, a spring barley experiment was monitored during three growing seasons in various campaigns with a TLS system and a field spectrometer. PH was derived from the TLS data and VIs from the hyperspectral data. Four major working tasks were carried out: (i) conduct extensive multi-annual field measurements during the growing seasons, (ii) derive bivariate biomass regression models (BRMs) from 3D spatial and spectral measurements for biomass estimations, (iii) fuse the 3D spatial and spectral data in multivariate BRMs to estimate biomass based on this extensive data set, and (iv) evaluate the robustness of the BRMs with a cross-validation.

2. Methods

2.1. Field Measurements

In three growing seasons (2012, 2013 and 2014), field experiments were carried out at the Campus Klein-Altendorf (50°37′51″N, E 6°59′32″) belonging to the Faculty of Agriculture at the University of Bonn, Germany. Due to crop rotation, the locations of the fields were slightly different between the years. However, soil and climatic conditions were similar with the surface of the soil being flat with a clayey silt luvisol and well-suited for crop cultivation [38]. According to the campus’ own weather records, the long-term average yearly precipitation was about 600 mm with a daily average temperature of 9.3 °C [39].

Each year, the field consisted of 36 small scale plots (3 × 7 m) where different cultivars of barley were cultivated with two levels of N fertilization. For half of the plots, a farmer’s common rate of 80 kg/ha∙N fertilizer was applied, for the other half a reduced rate of 40 kg/ha. In 2012 and 2013 each fertilization scheme was carried out once for 18 cultivars of spring barley (Barke, Wiebke, Beatrix, Eunova, Djamila, Streif, Ursa, Victoriana, Sissy, Perun, Apex, Isaria, Trumpf, Pflugs Intensiv, Heils Franken, Ackermanns Bavaria, Mauritia and Sebastian). In 2014, the set-up for the experiment was changed in that each fertilization scheme was repeated three times for six selected cultivars (Barke, Beatrix, Eunova, Trumpf, Mauritia and Sebastian). The experiments were carried out within the interdisciplinary research network CROP.SENSe.net (www.cropsense.uni-bonn.de). The research focus of this project was non-destructive sensor-based methods for detecting crop status such as nutrients, stress, and quality.

In this study, 3D spatial measurements from a TLS system, spectral measurements from a field spectrometer, and manual reference measurements were used. Due to the weather conditions the time of seeding changed and therefore so did the start of the growing season. The seeding dates were 21 March 2012, 9 April 2013, and 13 March 2014. In Table 1, all dates of TLS and spectrometer campaigns are listed as day after seeding (DAS) and a universal scale, known as the BBCH scale, was used to describe phenological stages and steps in the plant development, encoded in a decimal code [40,41]. The acronym BBCH is derived from the funding organizations: Biologische Bundesanstalt (German Federal Biological Research Centre for Agriculture and Forestry), Bundessortenamt (German Federal Office of Plant Varieties), and Chemical industry. The first number of the two-digit code represents the principal growth stage (Table 2) and the second subdivides further in short developmental steps. Through determining the BBCH codes during the growing seasons, the annual comparability was ensured. For each plot, the BBCH developmental step was determined as a mean of three plants. In Table 1, BBCH codes are given for the dates where plant parameters were manually measured. The codes are averaged per campaign, as the values were almost similar for all cultivars. Although the plant development varied among the years it can be seen that the BBCH codes indicate a comparable development.

As reference measurements, the heights of ten plants were measured for each plot and averaged in the post-processing. Moreover, in a defined sampling area of each plot, the above ground biomass of a 0.2 × 0.2 m area was destructively taken each time. The sampling area was neglected for the remote sensing measurements. In the laboratory, plants were cleaned and fresh weights were measured. After drying the samples for 120 h at 70 °C, dry biomass was weighted and extrapolated across the plot (g/m²).

2.1.1. Terrestrial Laser Scanning

The TLS configuration and setup was almost equal in all years. Thus for each campaign, the time-of-flight scanner Riegl LMS-Z420i was used (Figure 1A) [42]. The sensor operates with a near-infrared laser beam, has a beam divergence of 0.25 mrad, and a measurement rate of up to 11,000 points/sec. In addition its field of view is up to 80° in the vertical and 360° in the horizontal direction and this study used resolutions between 0.034° and 0.046°. The digital camera Nikon D200 was mounted on the laser scanner and the TLS point clouds were colorized from the images captured. Furthermore the sensor should be as high as possible above ground, resulting in a steep angle between scanner and investigated area enabling the best possible coverage of the crop surface and a homogenous penetration of the vegetation. Accordingly the scanner was mounted on the hydraulic platform of a tractor, raising the sensor to approximately 4 m above ground (Figure 1B). In order to lower shadowing effects and to attain an almost uniform spatial coverage, the field was scanned from its four corners. The coordinates of all scan positions and an additional target were required for the georeferencing and co-registration of the positions in the post-processing. Highly reflective cylinders arranged on ranging poles were used as targets (Figure 1C). These reflective cylinders can be easily detected by the scanner meaning their exact position in relation to the scan position can be measured [43]. The coordinates of the scan positions and ranging poles were measured with the highly accurate RTK-DGPS system Topcon HiPer Pro [44]. By establishing an own reference station each year, the precise merging of all data sets per year was ensured with the relative accuracy of this system being approximately 1 cm.

Table 1. Dates of the terrestrial laser scanning (TLS) and spectrometer (S) campaigns listed as day after seeding (DAS). Averaged codes for the developmental steps are given for the dates of manual plant parameter measurements (BBCH). For some dates BBCH codes were not determined (N/A).

	2012		2013		2014			2012		2013		2014			2012		2013		2014
DAS	Sensor	BBCH	Sensor	BBCH	Sensor	BBCH	DAS	Sensor	BBCH	Sensor	BBCH	Sensor	BBCH	DAS	Sensor	BBCH	Sensor	BBCH	Sensor	BBCH
15					TLS		45							75	TLS
16							46							76
17							47							77
18							48							78			TLS/ S	57
19							49			TLS/S	30			79
20			TLS/ S	N/A			50	TLS/ S						80
21							51							81
22							52							82					S
23							53							83		N/A
24							54		30			S		84	S				TLS	56
25							55							85
26							56					TLS	31	86	TLS
27							57							87
28							58	TLS						88
29							59							89
30							60							90
31							61							91			S	68
32							62							92			TLS
33							63							93
34			TLS				64			TLS/S	41			94
35			S	18			65							95
36							66							96						74
37							67							97					TLS/S
38							68							98
39							69							99
40							70	S	49			TLS/S	49	100
41					TLS/ S	29	71							101
42							72							102
43		N/A					73							103
44							74							104			TLS/ S	81

Table 1. Dates of the terrestrial laser scanning (TLS) and spectrometer (S) campaigns listed as day after seeding (DAS). Averaged codes for the developmental steps are given for the dates of manual plant parameter measurements (BBCH). For some dates BBCH codes were not determined (N/A).

	2012		2013		2014			2012		2013		2014			2012		2013		2014
DAS	Sensor	BBCH	Sensor	BBCH	Sensor	BBCH	DAS	Sensor	BBCH	Sensor	BBCH	Sensor	BBCH	DAS	Sensor	BBCH	Sensor	BBCH	Sensor	BBCH
15					TLS		45							75	TLS
16							46							76
17							47							77
18							48							78			TLS/ S	57
19							49			TLS/S	30			79
20			TLS/ S	N/A			50	TLS/ S						80
21							51							81
22							52							82					S
23							53							83		N/A
24							54		30			S		84	S				TLS	56
25							55							85
26							56					TLS	31	86	TLS
27							57							87
28							58	TLS						88
29							59							89
30							60							90
31							61							91			S	68
32							62							92			TLS
33							63							93
34			TLS				64			TLS/S	41			94
35			S	18			65							95
36							66							96						74
37							67							97					TLS/S
38							68							98
39							69							99
40							70	S	49			TLS/S	49	100
41					TLS/ S	29	71							101
42							72							102
43		N/A					73							103
44							74							104			TLS/ S	81

Furthermore, a digital terrain model (DTM) is required as a common reference surface for calculating plant height from the TLS data. In 2014, the bare ground of the field was scanned after seeding but before any vegetation was visible (Table 1: DAS 15). For technical reasons, it was not possible to acquire such data in 2012 and 2013, however, the ground was identifiable in the point cloud of the first campaigns due to the low and less dense vegetation.

Table 2. Principal growth stages of the BBCH scale.

Principal Growth Stage a	Stage Description	Time Frames Regarded for Biomass Estimation
0	Germination
1	Leaf development
2	Tillering		Pre-anthesis		Whole observed period
3	Stem elongation
4	Booting
5	Inflorescence emergence, heading
6	Flowering, anthesis
7	Development of fruit
8	Ripening
9	Senescence

^a first number of the two-digit code.

Table 2. Principal growth stages of the BBCH scale.

Principal Growth Stage a	Stage Description	Time Frames Regarded for Biomass Estimation
0	Germination
1	Leaf development
2	Tillering		Pre-anthesis		Whole observed period
3	Stem elongation
4	Booting
5	Inflorescence emergence, heading
6	Flowering, anthesis
7	Development of fruit
8	Ripening
9	Senescence

^a first number of the two-digit code.

Figure 1. Instrumental set-up: (A) terrestrial laser scanner Riegl LMS-Z420i; (B) tractor with hydraulic platform; (C) ranging pole with reflective cylinder.

Figure 1. Instrumental set-up: (A) terrestrial laser scanner Riegl LMS-Z420i; (B) tractor with hydraulic platform; (C) ranging pole with reflective cylinder.

2.1.2. Field Spectrometer Measurements

The ASD^® FieldSpec3 was used for measuring the reflectance several times during the growing seasons (all dates are listed in Table 1 above). This spectrometer measures the incoming light from 350 to 2500 nm with a sampling interval of 1.4 nm in the VNIR (350–1000 nm) and 2 nm in the SWIR (1001–2500 nm). These measurements are resampled to spectra with 1 nm resolution by the manufacturer’s software. At each position, ten measurements were taken and instantly averaged by the software, from 1 m above the canopy with a pistol grip, which was mounted on a cantilever to avoid shadows obscuring the sampling area. Additionally, a water level was used to ensure nadir view and no fore optic was used, resulting in a field of view of 25° and thus, a footprint area on the canopy with a radius of approximately 22 cm was achieved. Before the measurements, the spectrometer warmed up for at least 30 min and every 10 min or after illumination change, the spectrometer was optimized and calibrated with a spectralon calibration panel (polytetrafluoroethylene reference panel). Six positions were measured within each plot and for each position, the detector offset was corrected [16]. Then the six spectra were averaged, resulting in one spectrum per field plot, which was used in the further analysis.

2.2. Post-Processing

2.2.1. TLS Data

In the scanner software RiSCAN Pro, the DGPS data and the scans of all campaigns were imported into one project file per year. Based on the coordinates of the scan positions and reflectors, a direct georeferencing method was applied for the registration of all scan positions. However, a further adjustment was required due to small alignment errors between the point clouds. Based on the iterative closest point (ICP) algorithm [45], the Multi Station Adjustment in RiSCAN Pro allows the position and orientation of each scan position to be modified in multiple iterations and thus the best fitting result for all of them to be acquired. After optimizing the alignment with the ICP algorithm, the error, measured as standard deviation between used point-pairs, was 0.04 m on average for each campaign.

The point clouds were then merged to one dataset per campaign, and the area of interest was extracted. As reflections on insects or small particles in the air produced noise those points were manually removed. In addition a filtering scheme for selecting maximum points was used for determining the crop surface and in the same way, a filtering scheme for selecting minimum points was applied to extract ground points from the data sets of each first campaign. Finally, the data sets with XYZ coordinates of each point were exported.

The spatial analyses and visualization of the data were carried out in Esri ArcGIS Desktop 10.2.1. All point clouds were interpolated using the inverse distance weighting (IDW) algorithm, resulting in a raster with a consistent spatial resolution of 1 cm. IDW is an exact, deterministic algorithm that retains measured values at their sample location. The accuracy of measurements with a high density is maintained as all values are kept at their discrete location and not moved to fit the interpolation better [46]. As introduced by Hoffmeister et al. [43], the created raster data sets are referred to as crop surface models (CSMs). Similarly, a digital terrain model (DTM) was generated from the ground points and by subtracting the DTM from a CSM, plant heights were calculated pixel-wise. Moreover, by calculating the difference between two CSMs, plant growth was spatially measured. Hereinafter, growth is defined as temporal difference in height (for a detailed description of the CSMs creation and the calculation of plant heights see Tilly et al. [30]). The raster data sets with pixel-wise stored plant heights and growth were visualized as maps of plant height and growth, respectively. Then the plant heights were averaged plot-wise, allowing a common spatial base with the other measurements to be attained. It should be noted that previously, each plot was clipped with an inner buffer of 0.5 m to prevent border effects.

2.2.2. Spectral Data

For this study, established VIs were used to extract information from the hyperspectral data, measured with the field spectrometer. From the widespread of known hyper- and multispectral VIs for deriving different vegetation properties, six VIs were selected from the literature which have been found to be related to biomass and LAI. The selection was based on two criteria: Firstly, to make this study comparable to other studies VIs were selected which have been widely used in literature. Secondly, VIs with different spectral domains were used to examine if this would influence the prediction power of the fused models.

The NDVI was originally created for broad band satellite remote sensing [47] and has been widely used in the literature. It has been adapted to hyperspectral narrow bands and was specified for sensors such as GreenSeeker^TM and Crop Circle^TM [17]. Several articles reported relationships between the NDVI and biomass or LAI. However, NDVI has been shown to saturate in cases of dense and multi-layered canopy [13] and to have a non-linear relationship with biophysical parameters such as green LAI [48].

On this basis, Roujean and Breonin [49] developed the renormalized difference vegetation index (RDVI) for estimating the fraction of photosynthetically active radiation absorbed by vegetation, independent of a priori knowledge of the vegetation cover [49]. The RDVI showed strong relationships to LAI for different crops below an LAI of 5 [48,50]. In dense crop canopies with an LAI above five, RDVI tended to overestimate the LAI [48]. Simulations with the radiative transfer models PROSPECT and SAIL indicated that the RDVI is less affected by canopy structure, biochemistry, and soil background when estimating the LAI [50].

The red edge inflection point (REIP) was introduced by Guyot and Baret [51]. The REIP characterizes the inflection in the spectral red edge by calculating the wavelength with maximum slope. A variation of the inflection is mainly related to leaf chlorophyll content, leaf area index, and leaf inclination angle. Furthermore, soil reflectance and sun position have a limited effect [52].

GnyLi is a four-band VI for estimating biomass in the NIR and SWIR domain [17]. This VI was developed for winter wheat and showed good performance on different scales from plot to regional level and across several growth stages [17]. The GnyLi considers the two reflectance maxima and minima between 800 and 1300 nm. While the high reflectance is caused by the plants intercellular structure, the absorption at the minima is caused by cellulose, starch lignin, and water. These components contribute substantially to dry and fresh biomass and combining the two products helps to avoid saturation problems—this is a major advantage of this VI.

Similar to the GnyLi, the normalized reflectance index (NRI) was also developed for estimating biomass in winter wheat. The NRI was empirically developed by combining the shape of the NDVI and the best two band combination for biomass estimation with EO-1 Hyperion satellite data [53].

The red green blue vegetation index (RGBVI) was developed for estimating biomass based on bands available in a standard digital camera [54]. In this study, the RGB data was simulated from hyperspectral data where green, red, and blue values were calculated as the mean of the reflectance from 530 to 560 nm, 645 to 765 nm, and 465 to 495 nm, respectively. Thus, in contrast to other studies [37,54,55], the RGBVI was derived from radiometrically and spectrally calibrated data.

The six VIs used in this study can be categorized by the wavelength domains that are used in their formula. The NDVI, RDVI, and REIP use wavelengths in the visible and near-infrared domain (VISNIR VIs), the GnyLi and NRI use wavelengths in the near-infrared domain (NIR VIs), while the RGBVI uses wavelengths in the visible domain (VIS VI). The formulas of the VIs used in this study are given in Table 3, [47,49,51,53,54,56].

Table 3. Vegetation indices used in this study.

Wave-Length Domains	Vegetation Index	Formula	References
NIR	GnyLi	$\left(R900\times R1050-R955\times R1220\right)/(R900\times R1050+R955\times R1220)$	[56]
	NRI	$\left(R874-R1225\right)/\left(R874+R1225\right)$	[53]
VISNIR	NDVI	$\left(R798-R670\right)/\left(R798+R670\right)$	[47]
	RDVI	$\left(R798-R670\right)/\left(\sqrt{R798+R670}\right)$	[49]
	REIP	$700+40\ast \frac{\left(\frac{R670+R780}{2}\right)-R700}{R740-R700}$	[51]
VIS	RGBVI	$({R_{green}}^2-R_{blue}\times R_{red})/({R_{green}}^2+R_{blue}\times R_{red})$	[54]

Table 3. Vegetation indices used in this study.

Wave-Length Domains	Vegetation Index	Formula	References
NIR	GnyLi	$\left(R900\times R1050-R955\times R1220\right)/(R900\times R1050+R955\times R1220)$	[56]
	NRI	$\left(R874-R1225\right)/\left(R874+R1225\right)$	[53]
VISNIR	NDVI	$\left(R798-R670\right)/\left(R798+R670\right)$	[47]
	RDVI	$\left(R798-R670\right)/\left(\sqrt{R798+R670}\right)$	[49]
	REIP	$700+40\ast \frac{\left(\frac{R670+R780}{2}\right)-R700}{R740-R700}$	[51]
VIS	RGBVI	$({R_{green}}^2-R_{blue}\times R_{red})/({R_{green}}^2+R_{blue}\times R_{red})$	[54]

2.3. Biomass Regression Models

The main aim of this study was to establish biomass regression models (BRMs) and compare the potential of PH, VIs, and a fusion of PH and VIs for estimating barley biomass. The workflow for the BRM calibration and validation and the distinction of considered cases are shown in Figure 2. All calculations were performed in the R software environment [57]. The measurements from 2012 were excluded because the spectral data set was inconsistent, since due to unsuitable weather, no spectral data or only data for less than half of the plots could be acquired corresponding to the second and fourth TLS campaign, respectively (Table 1). Furthermore, as mentioned above, the number of cultivars was reduced in 2014 so as a result only these six cultivars were used from the 2013 data set to ensure comparability.

The reduced data set was split into four subsets to obtain independent values for calibration and validation. The first subset contained the plot-wise averaged measurements of plant height, calculated VIs and destructively taken biomass from 2013 (n = 48). Each other subset contained the same measurements of one repetition from 2014 (each n = 60). Thus, each subset contained the measurements of each cultivar from one plot with low and one with high N fertilizer level for the given campaign dates. A cross-validation was performed using these data sets: For each run, one subset was excluded from the BRM calibration and used for validating the resulting BRM.

First, bivariate BRMs for fresh and dry biomass were developed based on the CSM-derived PH or one of the six VIs. Linear and exponential BRMs were established since no trend regarding their usability for biomass estimations based on PH was clearly identifiable in earlier studies [31]. However, the biomass accumulation during the vegetative phase is exponential and other studies have shown that it is best estimated with exponential models [13,16]. For the exponential BRMs, the fresh and dry biomass values were natural log-transformed. Each BRM was calculated for two time frames, the whole observed period from tillering (BBCH stage 2) till the end of fruit development (BBCH stage 7) and the pre-anthesis period (till BBCH stage 6) (Table 2). The latter period is important as, for example, adequate crop conditions could buffer the grain yield against later environmental stress [10]. Thus, campaign numbers 3 to 6 and 2 to 6 were considered for 2013 and 2014, respectively, whereas each final campaign was excluded for the pre-anthesis BRMs. Considering the four possible subset combinations, overall 224 bivariate BRMs were established. Second, multivariate BRMs were established based on PH fused with each VI. Since they were also established as linear and exponential BRMs for fresh and dry biomass for both time frames, the four possible subset combinations led to 192 multivariate BRMs in total.

Figure 2. Workflow for the calibration and validation of the biomass regression models and distinction of cases for each model.

Figure 2. Workflow for the calibration and validation of the biomass regression models and distinction of cases for each model.

The calibration was evaluated by calculating the coefficient of determination (R²) for PH or VI vs. measured biomass and the standard error of the estimate (SE_E) [58]. For the validation, besides the R² (estimated vs. measured biomass), the root mean square error (RMSE), and Willmott’s index of agreement (d) [59,60] were determined. For each case, the results from the four runs were averaged. Finally, the robustness of the BRMs was evaluated by calculating the ratio between the R² values of BRM calibration and validation.

3. Results

3.1. Acquired Plant Parameters

The TLS-derived point clouds were used to establish CSMs and spatially calculate plant height. Results of the pixel-wise calculation were visualized in maps of plant height for each plot. As an example for this, maps of four plots and corresponding mean heights are shown in Figure 3 for the barley cultivar Trumpf. In the first campaign of 2013, plants were too small to obtain reasonable results. Thus, maps are presented for the last six and five campaigns of 2013 and 2014, respectively. One plot of each N fertilizer level is shown for both years. For the temporal development, an increase in plant height is observable until anthesis (BBCH stage 6) and afterwards, the development of ears begins and plant heights decrease due to the associated sinking of heads. Within all plots, the detailed representation of plant height is visible, which enables spatial differences in plant height to be detected. As a result, the exact calculation of mean heights can be assumed. A comparison of the plot-wise averaged values does not show that the fertilization rate directly influenced plant height.

Figure 3. Maps of four plots from the last six and five campaigns of 2013 and 2014, respectively. One plot of each N fertilizer level of the barley cultivar Trumpf is shown for each year (∅: Plot mean height).

Figure 3. Maps of four plots from the last six and five campaigns of 2013 and 2014, respectively. One plot of each N fertilizer level of the barley cultivar Trumpf is shown for each year (∅: Plot mean height).

The plot-wise averaged plant heights were used for statistical analysis and a comparison with the manual measurements. The linear regressions between all mean CSM-derived and manual measured plant heights for each of the three years is illustrated below in Figure 4. High coefficients of determination (R²) confirm the TLS-derived results. The R² across all years is 0.92, yearly separated values are also given in Figure 4. Moreover, a varying scattering between the years is indicated. The scattering is the lowest in the 2014 data set, which is presumably caused by the reduced number of cultivars in 2014 and associated with more similar plant heights.

Table A1. Statistics for the plot-wise averaged CSM-derived and manual measured plant heights (n: number of samples; $\overline{X}$: mean value; min: minimum; max: maximum; SD: standard deviation).

				CSM-Derived Plant Height (m)				Manual Measured Plant Height (m)
	BBCH	N level	n	$\overline{X}$	min	max	SD	$\overline{X}$	min	max	SD
2012
1	N/A	40	18	0.15	0.06	0.22	0.04	0.20	0.15	0.25	0.03
		80	18	0.18	0.14	0.24	0.03	0.20	0.15	0.26	0.03
2	30	40	18	0.21	0.13	0.28	0.04	0.35	0.28	0.42	0.05
		80	18	0.27	0.20	0.35	0.04	0.35	0.30	0.42	0.04
3	49	40	16	0.58	0.47	0.72	0.08	0.63	0.52	0.80	0.09
		80	15	0.64	0.48	0.80	0.11	0.66	0.54	0.79	0.08
4	N/A	40	14	0.73	0.61	0.81	0.06	0.86	0.74	0.96	0.06
		80	14	0.81	0.71	0.92	0.06	0.89	0.80	1.00	0.06
2013
2	18	40	18	0.21	0.05	0.37	0.09	0.19	0.14	0.25	0.03
		80	18	0.11	−0.07	0.25	0.08	0.20	0.16	0.27	0.03
3	30	40	18	0.33	0.15	0.51	0.11	0.29	0.19	0.56	0.09
		80	18	0.25	0.01	0.40	0.11	0.28	0.17	0.45	0.08
4	41	40	18	0.57	0.33	0.83	0.17	0.52	0.39	0.70	0.09
		80	18	0.56	0.24	0.79	0.18	0.57	0.31	0.81	0.13
5	57	40	16	0.84	0.64	1.11	0.13	0.77	0.66	0.95	0.07
		80	16	0.79	0.58	1.04	0.12	0.81	0.54	0.94	0.11
6	68	40	14	0.78	0.65	0.97	0.09	0.77	0.66	0.84	0.05
		80	14	0.77	0.66	0.90	0.08	0.83	0.76	1.00	0.06
7	81	40	14	0.75	0.62	0.96	0.10	0.72	0.65	0.82	0.06
		80	14	0.72	0.62	0.83	0.07	0.79	0.67	0.89	0.07
2014
2	29	40	18	0.16	0.12	0.24	0.03	0.19	0.12	0.30	0.04
		80	18	0.18	0.15	0.25	0.03	0.18	0.13	0.27	0.04
3	31	40	18	0.41	0.36	0.51	0.04	0.38	0.31	0.52	0.05
		80	18	0.42	0.34	0.52	0.05	0.36	0.27	0.45	0.05
4	49	40	18	0.63	0.53	0.70	0.04	0.59	0.53	0.65	0.03
		80	18	0.63	0.57	0.70	0.04	0.57	0.51	0.64	0.04
5	56	40	18	0.80	0.69	0.87	0.04	0.78	0.68	0.85	0.04
		80	18	0.81	0.75	0.93	0.04	0.78	0.72	0.89	0.04
6	74	40	18	0.76	0.66	0.84	0.04	0.77	0.68	0.83	0.03
		80	18	0.79	0.73	0.85	0.03	0.75	0.71	0.82	0.03

gives the mean, minimum, and maximum values of all plot-wise averaged values as well as the standard deviation per campaign of the CSM-derived and manual measured plant heights. Clearly observable lodging occurred in some plots between the second and third or fourth and fifth campaign in 2012 and 2013, respectively (for more details see Tilly et al. [61]). Those plots were neglected for the analysis and thus reduced the number of samples for the affected campaigns. As already stated for the visualized plots (Figure 3), an increase in plant height is detectable during pre-anthesis and a slight decrease is detectable afterwards. In addition, the difference between the mean values of both measurement methods is lower than 10% for almost all campaigns.

The field spectrometer measurements were used for calculating the six VIs (GnyLi, NDVI, NRI, RDVI, REIP, and RGBVI). As the spectral measurements from 2012 were not usable for a linkage with the TLS data, only the data sets from 2013 to 2014 were used. Moreover, from the data set of 2013 only measurements of the cultivars selected in 2014 were considered and the data sets of plant height and biomass were accordingly adapted to ensure comparability. For each campaign, the values for both N fertilizer levels were averaged. Table 4 shows the statistics for the reduced data sets of the nine regarded campaigns. Additionally, the yearly mean biomass values were calculated for the pre-anthesis and whole observed period, as reference for the later evaluation of the biomass estimation.

Table 4. Statistics for the plot-wise averaged CSM-derived plant heights and destructively taken biomass for the reduced data sets of 2013 and 2014 (n: number of samples; $\overline{X}$: mean value; min: minimum; max: maximum; SD: standard deviation).

		CSM-Derived Plant Height (m)				Fresh Biomass (g/m2)				Dry Biomass (g/m2)
	n	$\overline{X}$	min	max	SD	$\overline{X}$	min	max	SD	$\overline{X}$	min	max	SD
2013
3	12	0.22	0.01	0.39	0.13	51.32	19.64	86.90	19.77	6.73	2.08	10.88	2.36
4	12	0.47	0.24	0.71	0.17	115.66	62.41	178.62	33.68	16.63	8.20	29.00	6.09
5	12	0.78	0.58	0.99	0.13	202.82	106.75	309.20	65.24	35.34	17.38	57.17	13.74
6	12	0.78	0.65	0.93	0.07	185.27	119.45	306.23	49.88	50.36	35.47	67.50	9.19
	Mean pre-anthesis period					123.27	62.93	191.57	39.56	19.57	9.22	32.35	7.40
	Mean whole observed period					138.77	77.06	220.24	42.14	27.27	15.78	41.14	7.85
2014
2	36	0.17	0.12	0.25	0.03	26.25	10.65	44.66	8.20	3.56	1.32	6.21	1.12
3	36	0.41	0.34	0.52	0.04	89.08	49.07	129.46	21.57	11.59	6.63	16.71	2.68
4	36	0.63	0.53	0.70	0.04	113.02	65.75	166.48	24.47	18.62	11.06	28.27	3.97
5	36	0.81	0.69	0.99	0.05	127.41	84.26	217.33	27.90	31.09	19.45	50.85	6.33
6	36	0.78	0.66	0.99	0.05	142.77	79.79	241.76	36.45	46.66	26.10	75.05	11.21
	Mean pre-anthesis period					88.94	52.43	139.48	20.54	16.22	9.62	25.51	3.53
	Mean whole observed period					99.71	57.90	159.94	23.72	22.30	12.91	35.42	5.06

Table 4. Statistics for the plot-wise averaged CSM-derived plant heights and destructively taken biomass for the reduced data sets of 2013 and 2014 (n: number of samples; $\overline{X}$: mean value; min: minimum; max: maximum; SD: standard deviation).

		CSM-Derived Plant Height (m)				Fresh Biomass (g/m2)				Dry Biomass (g/m2)
	n	$\overline{X}$	min	max	SD	$\overline{X}$	min	max	SD	$\overline{X}$	min	max	SD
2013
3	12	0.22	0.01	0.39	0.13	51.32	19.64	86.90	19.77	6.73	2.08	10.88	2.36
4	12	0.47	0.24	0.71	0.17	115.66	62.41	178.62	33.68	16.63	8.20	29.00	6.09
5	12	0.78	0.58	0.99	0.13	202.82	106.75	309.20	65.24	35.34	17.38	57.17	13.74
6	12	0.78	0.65	0.93	0.07	185.27	119.45	306.23	49.88	50.36	35.47	67.50	9.19
	Mean pre-anthesis period					123.27	62.93	191.57	39.56	19.57	9.22	32.35	7.40
	Mean whole observed period					138.77	77.06	220.24	42.14	27.27	15.78	41.14	7.85
2014
2	36	0.17	0.12	0.25	0.03	26.25	10.65	44.66	8.20	3.56	1.32	6.21	1.12
3	36	0.41	0.34	0.52	0.04	89.08	49.07	129.46	21.57	11.59	6.63	16.71	2.68
4	36	0.63	0.53	0.70	0.04	113.02	65.75	166.48	24.47	18.62	11.06	28.27	3.97
5	36	0.81	0.69	0.99	0.05	127.41	84.26	217.33	27.90	31.09	19.45	50.85	6.33
6	36	0.78	0.66	0.99	0.05	142.77	79.79	241.76	36.45	46.66	26.10	75.05	11.21
	Mean pre-anthesis period					88.94	52.43	139.48	20.54	16.22	9.62	25.51	3.53
	Mean whole observed period					99.71	57.90	159.94	23.72	22.30	12.91	35.42	5.06

Figure 4. Regression of the mean CSM-derived and manual measured plant heights (2012: n = 131; 2013: n = 196; 2014: n = 180.).

Figure 4. Regression of the mean CSM-derived and manual measured plant heights (2012: n = 131; 2013: n = 196; 2014: n = 180.).

3.2. Biomass Estimation

The barley biomass was estimated by establishing 224 bivariate and 192 multivariate biomass regression models (BRMs) based on plant height (PH) and vegetation indices (VIs). Table 5 shows the statistical parameters for the BRM calibration. The table is vertically divided into bivariate or multivariate BRMs and the regarded time frames. Horizontally it distinguishes between dry or fresh biomass and linear or exponential BRMs. However, the results of the linear and exponential BRMs cannot be directly compared due to the log-transformation of biomass for the latter ones. Since the biomass accumulation during the vegetative phase is exponential and other studies have shown that it is best estimated with exponential BRMs [13,16] only the exponential BRMs are regarded in the following. For each model the coefficient of determination (R²) and the standard error of the estimate (SE_E) are given as mean values of the four possible subset combinations.

Each established BRM was validated with the remaining fourth subsets. Table 6 shows the R², root mean square error (RMSE), and Willmott’s index of agreement (d) for the model validation as mean values of the four subset combinations. The subdivision of the table is equivalent to that of Table 5. The results of the bivariate BRMs are regarded in the following subsection; the fusion of both plant parameters to multivariate BRMs is examined in the last subsection of this chapter. As the results of the calibration and validation show a similar tendency, only the values of the validation are stated. However, to evaluate the robustness of the BRMs, an overall comparison of differences between calibration and validation is given at the end of this chapter.

Table 5. Statistics for the model calibration as mean values of the four subset combinations (R²: coefficient of determination; SE_E: standard error of the estimate).

		Bivariate BRMs					Multivariate BRMs
			Whole Period		Pre-Anthesis			Whole Period		Pre-Anthesis
		Estimator	R2	${{SE}_E}^a$	R2	${{SE}_E}^a$	Estimatorb	R2	${{SE}_E}^a$	R2	${{SE}_E}^a$
Dry biomass	Linear	PH	0.65	10.03	0.76	5.73
		GnyLi	0.52	11.75	0.68	6.67	GnyLi	0.65	34.63	0.77	25.41
		NDVI	0.07	16.38	0.34	9.56	NDVI	0.69	21.49	0.76	20.73
		NRI	0.54	11.58	0.70	6.40	NRI	0.65	35.04	0.77	24.86
		RDVI	0.13	15.87	0.39	9.21	RDVI	0.69	19.18	0.76	21.40
		REIP	0.12	15.92	0.58	7.60	REIP	0.73	1933.86	0.76	258.29
		RGBVI	0.05	16.55	0.26	10.10	RGBVI	0.68	22.28	0.76	23.23
	Exponential	PH	0.84	0.37	0.84	0.34
		GnyLi	0.80	0.42	0.85	0.32	GnyLi	0.86	2.43	0.88	2.14
		NDVI	0.30	0.77	0.61	0.53	NDVI	0.85	2.84	0.88	3.99
		NRI	0.81	0.40	0.87	0.30	NRI	0.87	2.29	0.89	1.96
		RDVI	0.41	0.71	0.68	0.48	RDVI	0.85	2.52	0.88	2.84
		REIP	0.37	0.73	0.77	0.40	REIP	0.84	30.37	0.86	48.49
		RGBVI	0.23	0.81	0.48	0.60	RGBVI	0.85	2.51	0.87	2.73
		Estimator	R2	${{\mathit{S}\mathit{E}}_{\mathit{E}}}^{\mathit{a}}$	R2	${{\mathit{S}\mathit{E}}_{\mathit{E}}}^{\mathit{a}}$	Estimator b	R2	${{\mathit{S}\mathit{E}}_{\mathit{E}}}^{\mathit{a}}$	R2	${{\mathit{S}\mathit{E}}_{\mathit{E}}}^{\mathit{a}}$
Fresh biomass	Linear	PH	0.59	36.08	0.60	33.73
		GnyLi	0.58	36.55	0.62	33.18	GnyLi	0.62	131.81	0.64	118.76
		NDVI	0.25	48.90	0.42	40.91	NDVI	0.60	182.47	0.63	200.34
		NRI	0.59	36.40	0.62	32.85	NRI	0.62	122.25	0.64	108.72
		RDVI	0.35	45.74	0.50	37.81	RDVI	0.61	152.56	0.64	158.23
		REIP	0.30	47.23	0.55	35.78	REIP	0.60	584.00	0.63	2366.77
		RGBVI	0.22	49.75	0.37	42.66	RGBVI	0.61	160.32	0.64	155.26
	Exponential	PH	0.70	0.37	0.68	0.39
		GnyLi	0.76	0.33	0.76	0.34	GnyLi	0.77	1.86	0.77	1.77
		NDVI	0.46	0.50	0.65	0.41	NDVI	0.77	3.74	0.79	4.30
		NRI	0.77	0.33	0.77	0.33	NRI	0.77	1.66	0.77	1.55
		RDVI	0.59	0.43	0.74	0.35	RDVI	0.79	2.69	0.82	2.89
		REIP	0.47	0.49	0.71	0.37	REIP	0.72	22.30	0.74	73.14
		RGBVI	0.38	0.53	0.55	0.46	RGBVI	0.77	2.58	0.78	2.68

^a The SE_E for exponential models is calculated from natural log-transformed biomass values; ^b each fused with PH.

Table 5. Statistics for the model calibration as mean values of the four subset combinations (R²: coefficient of determination; SE_E: standard error of the estimate).

		Bivariate BRMs					Multivariate BRMs
			Whole Period		Pre-Anthesis			Whole Period		Pre-Anthesis
		Estimator	R2	${{SE}_E}^a$	R2	${{SE}_E}^a$	Estimatorb	R2	${{SE}_E}^a$	R2	${{SE}_E}^a$
Dry biomass	Linear	PH	0.65	10.03	0.76	5.73
		GnyLi	0.52	11.75	0.68	6.67	GnyLi	0.65	34.63	0.77	25.41
		NDVI	0.07	16.38	0.34	9.56	NDVI	0.69	21.49	0.76	20.73
		NRI	0.54	11.58	0.70	6.40	NRI	0.65	35.04	0.77	24.86
		RDVI	0.13	15.87	0.39	9.21	RDVI	0.69	19.18	0.76	21.40
		REIP	0.12	15.92	0.58	7.60	REIP	0.73	1933.86	0.76	258.29
		RGBVI	0.05	16.55	0.26	10.10	RGBVI	0.68	22.28	0.76	23.23
	Exponential	PH	0.84	0.37	0.84	0.34
		GnyLi	0.80	0.42	0.85	0.32	GnyLi	0.86	2.43	0.88	2.14
		NDVI	0.30	0.77	0.61	0.53	NDVI	0.85	2.84	0.88	3.99
		NRI	0.81	0.40	0.87	0.30	NRI	0.87	2.29	0.89	1.96
		RDVI	0.41	0.71	0.68	0.48	RDVI	0.85	2.52	0.88	2.84
		REIP	0.37	0.73	0.77	0.40	REIP	0.84	30.37	0.86	48.49
		RGBVI	0.23	0.81	0.48	0.60	RGBVI	0.85	2.51	0.87	2.73
		Estimator	R2	${{\mathit{S}\mathit{E}}_{\mathit{E}}}^{\mathit{a}}$	R2	${{\mathit{S}\mathit{E}}_{\mathit{E}}}^{\mathit{a}}$	Estimator b	R2	${{\mathit{S}\mathit{E}}_{\mathit{E}}}^{\mathit{a}}$	R2	${{\mathit{S}\mathit{E}}_{\mathit{E}}}^{\mathit{a}}$
Fresh biomass	Linear	PH	0.59	36.08	0.60	33.73
		GnyLi	0.58	36.55	0.62	33.18	GnyLi	0.62	131.81	0.64	118.76
		NDVI	0.25	48.90	0.42	40.91	NDVI	0.60	182.47	0.63	200.34
		NRI	0.59	36.40	0.62	32.85	NRI	0.62	122.25	0.64	108.72
		RDVI	0.35	45.74	0.50	37.81	RDVI	0.61	152.56	0.64	158.23
		REIP	0.30	47.23	0.55	35.78	REIP	0.60	584.00	0.63	2366.77
		RGBVI	0.22	49.75	0.37	42.66	RGBVI	0.61	160.32	0.64	155.26
	Exponential	PH	0.70	0.37	0.68	0.39
		GnyLi	0.76	0.33	0.76	0.34	GnyLi	0.77	1.86	0.77	1.77
		NDVI	0.46	0.50	0.65	0.41	NDVI	0.77	3.74	0.79	4.30
		NRI	0.77	0.33	0.77	0.33	NRI	0.77	1.66	0.77	1.55
		RDVI	0.59	0.43	0.74	0.35	RDVI	0.79	2.69	0.82	2.89
		REIP	0.47	0.49	0.71	0.37	REIP	0.72	22.30	0.74	73.14
		RGBVI	0.38	0.53	0.55	0.46	RGBVI	0.77	2.58	0.78	2.68

^a The SE_E for exponential models is calculated from natural log-transformed biomass values; ^b each fused with PH.

3.2.1. Bivariate Models

All cases show moderate to good results for bivariate BRMs based on PH. For each time frame, PH shows the same and similar relationship with dry and fresh biomass, respectively (Table 6). Scatterplots of measured vs. estimated biomass for selected examples are shown in the last subsection in comparison with multivariate BRMs.

Table 6. Statistics for the model validation as mean values of the four subset combinations (R²: coefficient of determination; RMSE: root mean square error (g/m²); d: Willmott’s index of agreement).

		Bivariate BRMs							Multivariate BRMs
			Whole Period			Pre-Anthesis				Whole Period			Pre-Anthesis
		Estimator	R2	RMSE a	d	R2	RMSE a	d	Estimator b	R2	RMSE a	d	R2	RMSE a	d
Dry biomass	Linear	PH	0.66	10.30	0.88	0.80	5.91	0.92
		GnyLi	0.54	11.99	0.81	0.72	6.93	0.88	GnyLi	0.65	10.49	0.88	0.79	5.93	0.92
		NDVI	0.07	16.51	0.33	0.38	9.78	0.64	NDVI	0.71	10.01	0.89	0.80	5.93	0.92
		NRI	0.55	11.82	0.82	0.74	6.66	0.89	NRI	0.66	10.47	0.88	0.80	5.91	0.92
		RDVI	0.13	16.01	0.44	0.41	9.34	0.71	RDVI	0.72	9.89	0.89	0.80	5.93	0.92
		REIP	0.15	16.20	0.46	0.68	7.90	0.83	REIP	0.73	9.14	0.91	0.80	5.91	0.92
		RGBVI	0.04	16.66	0.26	0.28	10.18	0.58	RGBVI	0.70	10.45	0.88	0.80	5.97	0.92
	Exponential	PH	0.85	0.39	0.95	0.85	0.36	0.94
		GnyLi	0.80	0.42	0.94	0.86	0.33	0.95	GnyLi	0.87	0.36	0.96	0.89	0.31	0.96
		NDVI	0.29	0.77	0.63	0.59	0.54	0.81	NDVI	0.85	0.38	0.95	0.87	0.30	0.96
		NRI	0.81	0.40	0.94	0.87	0.31	0.96	NRI	0.87	0.36	0.96	0.89	0.29	0.96
		RDVI	0.40	0.71	0.73	0.66	0.48	0.87	RDVI	0.85	0.38	0.95	0.88	0.30	0.96
		REIP	0.40	0.75	0.72	0.82	0.43	0.90	REIP	0.85	0.39	0.95	0.89	0.34	0.95
		RGBVI	0.22	0.82	0.55	0.48	0.62	0.75	RGBVI	0.85	0.38	0.95	0.86	0.31	0.96
		Estimator	R2	RMSEa	d	R2	RMSEa	d	Estimator b	R2	RMSEa	d	R2	RMSEa	d
Fresh biomass	Linear	PH	0.67	38.54	0.84	0.70	35.70	0.85
		GnyLi	0.65	38.83	0.83	0.72	35.45	0.84	GnyLi	0.69	37.59	0.85	0.74	34.47	0.86
		NDVI	0.27	50.16	0.58	0.51	42.15	0.70	NDVI	0.67	38.10	0.84	0.73	34.51	0.85
		NRI	0.65	38.50	0.83	0.72	34.95	0.85	NRI	0.69	37.54	0.85	0.74	34.32	0.86
		RDVI	0.38	47.01	0.67	0.59	38.58	0.77	RDVI	0.68	37.76	0.85	0.74	33.65	0.86
		REIP	0.41	49.76	0.66	0.77	38.07	0.81	REIP	0.67	38.67	0.84	0.77	36.35	0.84
		RGBVI	0.21	50.41	0.53	0.41	42.65	0.67	RGBVI	0.66	37.96	0.85	0.71	34.12	0.86
	Exponential	PH	0.73	0.40	0.89	0.71	0.42	0.88
		GnyLi	0.77	0.35	0.92	0.79	0.36	0.91	GnyLi	0.79	0.34	0.92	0.80	0.36	0.92
		NDVI	0.44	0.51	0.73	0.64	0.42	0.83	NDVI	0.78	0.34	0.92	0.79	0.34	0.92
		NRI	0.77	0.34	0.92	0.79	0.35	0.92	NRI	0.79	0.34	0.92	0.79	0.35	0.92
		RDVI	0.57	0.44	0.82	0.73	0.36	0.89	RDVI	0.80	0.33	0.93	0.83	0.31	0.93
		REIP	0.54	0.53	0.77	0.82	0.42	0.87	REIP	0.77	0.39	0.90	0.82	0.40	0.88
		RGBVI	0.36	0.54	0.68	0.53	0.47	0.78	RGBVI	0.76	0.34	0.92	0.76	0.34	0.92

^a The RMSE for exponential models is calculated from natural log-transformed biomass values; ^b each fused with PH.

Table 6. Statistics for the model validation as mean values of the four subset combinations (R²: coefficient of determination; RMSE: root mean square error (g/m²); d: Willmott’s index of agreement).

		Bivariate BRMs							Multivariate BRMs
			Whole Period			Pre-Anthesis				Whole Period			Pre-Anthesis
		Estimator	R2	RMSE a	d	R2	RMSE a	d	Estimator b	R2	RMSE a	d	R2	RMSE a	d
Dry biomass	Linear	PH	0.66	10.30	0.88	0.80	5.91	0.92
		GnyLi	0.54	11.99	0.81	0.72	6.93	0.88	GnyLi	0.65	10.49	0.88	0.79	5.93	0.92
		NDVI	0.07	16.51	0.33	0.38	9.78	0.64	NDVI	0.71	10.01	0.89	0.80	5.93	0.92
		NRI	0.55	11.82	0.82	0.74	6.66	0.89	NRI	0.66	10.47	0.88	0.80	5.91	0.92
		RDVI	0.13	16.01	0.44	0.41	9.34	0.71	RDVI	0.72	9.89	0.89	0.80	5.93	0.92
		REIP	0.15	16.20	0.46	0.68	7.90	0.83	REIP	0.73	9.14	0.91	0.80	5.91	0.92
		RGBVI	0.04	16.66	0.26	0.28	10.18	0.58	RGBVI	0.70	10.45	0.88	0.80	5.97	0.92
	Exponential	PH	0.85	0.39	0.95	0.85	0.36	0.94
		GnyLi	0.80	0.42	0.94	0.86	0.33	0.95	GnyLi	0.87	0.36	0.96	0.89	0.31	0.96
		NDVI	0.29	0.77	0.63	0.59	0.54	0.81	NDVI	0.85	0.38	0.95	0.87	0.30	0.96
		NRI	0.81	0.40	0.94	0.87	0.31	0.96	NRI	0.87	0.36	0.96	0.89	0.29	0.96
		RDVI	0.40	0.71	0.73	0.66	0.48	0.87	RDVI	0.85	0.38	0.95	0.88	0.30	0.96
		REIP	0.40	0.75	0.72	0.82	0.43	0.90	REIP	0.85	0.39	0.95	0.89	0.34	0.95
		RGBVI	0.22	0.82	0.55	0.48	0.62	0.75	RGBVI	0.85	0.38	0.95	0.86	0.31	0.96
		Estimator	R2	RMSEa	d	R2	RMSEa	d	Estimator b	R2	RMSEa	d	R2	RMSEa	d
Fresh biomass	Linear	PH	0.67	38.54	0.84	0.70	35.70	0.85
		GnyLi	0.65	38.83	0.83	0.72	35.45	0.84	GnyLi	0.69	37.59	0.85	0.74	34.47	0.86
		NDVI	0.27	50.16	0.58	0.51	42.15	0.70	NDVI	0.67	38.10	0.84	0.73	34.51	0.85
		NRI	0.65	38.50	0.83	0.72	34.95	0.85	NRI	0.69	37.54	0.85	0.74	34.32	0.86
		RDVI	0.38	47.01	0.67	0.59	38.58	0.77	RDVI	0.68	37.76	0.85	0.74	33.65	0.86
		REIP	0.41	49.76	0.66	0.77	38.07	0.81	REIP	0.67	38.67	0.84	0.77	36.35	0.84
		RGBVI	0.21	50.41	0.53	0.41	42.65	0.67	RGBVI	0.66	37.96	0.85	0.71	34.12	0.86
	Exponential	PH	0.73	0.40	0.89	0.71	0.42	0.88
		GnyLi	0.77	0.35	0.92	0.79	0.36	0.91	GnyLi	0.79	0.34	0.92	0.80	0.36	0.92
		NDVI	0.44	0.51	0.73	0.64	0.42	0.83	NDVI	0.78	0.34	0.92	0.79	0.34	0.92
		NRI	0.77	0.34	0.92	0.79	0.35	0.92	NRI	0.79	0.34	0.92	0.79	0.35	0.92
		RDVI	0.57	0.44	0.82	0.73	0.36	0.89	RDVI	0.80	0.33	0.93	0.83	0.31	0.93
		REIP	0.54	0.53	0.77	0.82	0.42	0.87	REIP	0.77	0.39	0.90	0.82	0.40	0.88
		RGBVI	0.36	0.54	0.68	0.53	0.47	0.78	RGBVI	0.76	0.34	0.92	0.76	0.34	0.92

^a The RMSE for exponential models is calculated from natural log-transformed biomass values; ^b each fused with PH.

Most VIs lead to better results for pre-anthesis than for the whole observed period. For dry biomass, the RGBVI performs worst for both time frames (Table 6, top left quarter). The largest difference between the whole observed period and the pre-anthesis can be found for the NDVI (R² = 0.29 vs. 0.59), while the NIR VIs as the GnyLi perform more consistently (R² = 0.80 vs. 0.86). Both, the NRI and the GnyLi also reveal best results for pre-anthesis (R² = 0.87, 0.86) and for the whole observed period (R² = 0.81, 0.80). In pre-anthesis, the relative difference between the NIR VIs and VISNIR VIs is smaller. Figure 5 shows scatterplots of measured vs. estimated dry biomass of one validation dataset for selected VIs and as expected from the high R² values, the estimated biomass from the GnyLi BRM corresponds well with the measured biomass (close to the 1:1 line). In pre-anthesis, the same applies the REIP whereas the NDVI and RGBVI saturate at about 7.4 g/m². For the whole observed period, biomass estimated by the BRM of REIP, NDVI and RGBVI does not align well with what was measured. The scatterplots reveal that the dynamic range of the models does not cover the range of the measured biomass values.

Figure 5. Scatterplots of measured vs. estimated dry biomass for one validation data set for NDVI, RGBVI REIP, and GnyLi (exponential model). Pre-anthesis: crosses and solid green line; whole observed period: circles and dashed black line; 1:1 line: light grey.

Figure 5. Scatterplots of measured vs. estimated dry biomass for one validation data set for NDVI, RGBVI REIP, and GnyLi (exponential model). Pre-anthesis: crosses and solid green line; whole observed period: circles and dashed black line; 1:1 line: light grey.

Better results are also obtained for pre-anthesis of fresh biomass than for the whole observed period, although the differences are smaller than for dry biomass. The NIR VIs perform most consistently for both periods and have the highest R² values for the whole observed period. However, particularly for the whole observed period, the relative difference between the NIR VIs and the VIS and VISNIR VIs is smaller than for dry biomass and in pre-anthesis, the relative difference between the NIR VIs and other VIs is further reduced. Additionally, the REIP (R² = 0.82) yields better results than the NIR VIs (each R² = 0.79). Again, the RGBVI performs worst. Figure 6 shows scatterplots of measured vs. estimated fresh biomass of one validation dataset for selected VIs. As expected from the high R² values, biomass estimated from the GnyLi BRM corresponds well with the measured values (close to the 1:1 line). In pre-anthesis, the same applies for the REIP, whereas the NDVI and RGBVI saturate at about 55 g/m². As for dry biomass, the BRMs based on the REIP and particularly the NDVI and RGBVI show a poor relationship between estimated and measured fresh biomass. Overall, most VISNIR VIs and the RGBVI yield better results for fresh biomass than for dry biomass. The NIR VIs perform best and most consistently (Table 6, bottom left quarter).

Figure 6. Scatterplots of measured vs. estimated fresh biomass for one validation data set for NDVI, RGBVI REIP, and GnyLi (exponential model). Pre-anthesis: crosses and solid green line; whole observed period: circles and dashed black line; 1:1 line: light grey.

Figure 6. Scatterplots of measured vs. estimated fresh biomass for one validation data set for NDVI, RGBVI REIP, and GnyLi (exponential model). Pre-anthesis: crosses and solid green line; whole observed period: circles and dashed black line; 1:1 line: light grey.

3.2.2. Multivariate Models

For dry biomass, PH is the best individual estimator across the whole observed period (R² = 0.85) and a slight improvement is only achieved when fused with one of the NIR VIs in a multivariate BRM (both R² = 0.87). In pre-anthesis, PH and the NIR VIs perform similarly to the bivariate BRMs (R² = 0.85, 0.86, 0.87) and when PH is fused with the NIR VIs or the REIP, the predictability slightly increases (R² = 0.89).

For fresh biomass across the whole observed period, PH (R² = 0.73) yields comparable results to the NIR VIs (both R² = 0.77) although the fusion of PH with NIR VIs slightly improves the estimation (both R² = 0.79). Only the multivariate BRM from PH and RDVI is very slightly better (R² = 0.80). In pre-anthesis, REIP, GnyLi, NRI, and RDVI explain up to 11% more variation (R² = 0.82, 0.79, 0.79, 0.73) then PH (R² = 0.71). When PH is fused with any VI, the predictability is improved compared to most individual estimators and even the RGBVI in combination with PH improves the estimation of dry and fresh biomass for pre-anthesis yielding an R² of 0.71 and 0.76, respectively. In the fused analysis, the RGBVI performs only slightly weaker than the other VIs. Nevertheless, only the RDVI fused with PH slightly increases the predictability (R² = 0.83) compared to the bivariate BRM based on the RDVI. Figure 7 shows the scatterplots of measured vs. estimated values of one validation dataset from the bivariate BRM of PH and the multivariate BRM of PH and GnyLi for dry biomass in pre-anthesis and fresh biomass across the whole observed period. The model fit is only slightly improved by fusing PH with the VI.

Figure 7. Scatterplot for one validation data set for the pre-anthesis (green) and for the whole observed period (black) of the bivariate BRM of PH (circles and solid regression line) and multivariate BRM of PH and GnyLi (crosses and dashed regression line) for fresh biomass (top) and dry biomass (bottom) (all exponential models); 1:1 line: light grey.

Figure 7. Scatterplot for one validation data set for the pre-anthesis (green) and for the whole observed period (black) of the bivariate BRM of PH (circles and solid regression line) and multivariate BRM of PH and GnyLi (crosses and dashed regression line) for fresh biomass (top) and dry biomass (bottom) (all exponential models); 1:1 line: light grey.

The robustness of the models was evaluated by calculating the ratio between the R² values of model calibration and validation for each BRM (

Table A2. Ratio between model calibration and validation ($R_{cal}^2$: coefficient of determination from calibration; $R_{val}^2$: coefficient of determination from validation).

		Bivariate BRMs			Multivariate BRMs
			Whole Period	Pre-Anthesis		Whole Period	Pre-Anthesis
		Estimator	$R_{cal}^2$/$R_{val}^2$	$R_{cal}^2$/$R_{val}^2$	Estimatora	$R_{cal}^2$/$R_{val}^2$	$R_{cal}^2$/$R_{val}^2$
Dry biomass	Linear	PH	0.98	0.95
		GnyLi	0.96	0.94	GnyLi	1.00	0.97
		NDVI	1.00	0.89	NDVI	0.97	0.95
		NRI	0.98	0.95	NRI	0.98	0.96
		RDVI	1.00	0.95	RDVI	0.96	0.95
		REIP	0.80	0.85	REIP	1.00	0.95
		RGBVI	1.25	0.93	RGBVI	0.97	0.95
	Exponential	PH	0.99	0.99
		GnyLi	1.00	0.99	GnyLi	0.99	0.99
		NDVI	1.03	1.03	NDVI	1.00	1.01
		NRI	1.00	1.00	NRI	1.00	1.00
		RDVI	1.03	1.03	RDVI	1.00	1.00
		REIP	0.93	0.94	REIP	0.99	0.97
		RGBVI	1.05	1.00	RGBVI	1.00	1.01
Fresh biomass	Linear	PH	0.88	0.86
		GnyLi	0.89	0.86	GnyLi	0.90	0.86
		NDVI	0.93	0.82	NDVI	0.90	0.86
		NRI	0.91	0.86	NRI	0.90	0.86
		RDVI	0.92	0.85	RDVI	0.90	0.86
		REIP	0.73	0.71	REIP	0.90	0.82
		RGBVI	1.05	0.90	RGBVI	0.92	0.90
	Exponential	PH	0.96	0.96
		GnyLi	0.99	0.96	GnyLi	0.97	0.96
		NDVI	1.05	1.02	NDVI	0.99	1.00
		NRI	1.00	0.97	NRI	0.97	0.97
		RDVI	1.04	1.01	RDVI	0.99	0.99
		REIP	0.87	0.87	REIP	0.94	0.90
		RGBVI	1.06	1.04	RGBVI	1.01	1.03

^a each fused with PH.

). Since the R² of calibration was divided through the R² of validation, values above 1 indicate better results from the calibration and below 1 indicate better results from the validation. Consequently, values close to 1 show a robust performance. For the bivariate BRMs, PH and almost all VIs are supported as robust estimators by ratios close to 1 for all cases. The weakest ratios are attained for the REIP, in particular for fresh biomass with linear BRMs (0.73, 0.71). For the multivariate BRMs, good ratios are found for all cases. Only the linear BRMs for fresh biomass show slightly weaker values for the pre-anthesis period.

4. Discussion

The overall aim of this study was to evaluate whether the fusion of PH and VIs can improve the predictability of dry and fresh barley biomass compared to each parameter as individual estimator. For this work, the use of TLS to derive PH was verified and bivariate BRMs based on PH or one of six VIs as well as multivariate BRMs based on the fusion of PH with each VI were established. Extensive fieldwork over three years supported the practical application of the presented methods for monitoring crop development on plot level. The same instruments were used for all measurements whereby variations through different sensors could be excluded. However, the design of the field experiment and the measurement program was slightly modified and optimized over the years. Hence, only a part of the acquired data was used for the final model generation in order to ensure the comparability between the data sets. In the following, first the retrieval of PH from TLS data is discussed before the different BRMs are examined.

4.1. TLS-Derived Plant Height

The presented study verified the reliability of the laser scanner Riegl LMS-Z420i for capturing crop surfaces. In comparison with past studies [31,43], the scanning angle to the field was optimized through the elevated position on the hydraulic platform. However, uncertainties still remain about the influence of the scanning angle and the fixed position of the scanner during the measurements. As maintained by Ehlert and Heisig [62]—the scanning angle can cause overestimations in the height of reflection points and should be considered in the calculation of heights. In this study, the crop surface was determined from the merged and cleaned point clouds of four scan positions, filtered with a scheme for selecting maximum points. Overestimations should therefore be precluded.

For the practical implementation of CSM-derived plant height measurements, further aspects have to be considered. Usually, the factors time and cost have a major influence on choosing a system. As shown by Hämmerle and Höfle [63] the appropriate point density for generating a CSM varies depending on the application. In further studies, cost-efficient systems, such as the Velodyne HDL-64E [64], should be considered to investigate their potential for capturing crop surfaces in an adequate resolution. In the distant future, low-cost stationary systems might get permanently established for monitoring plant growth on field level. Moreover, recent developments have brought up new laser scanning platforms that might accelerate the field measurement process and optimize the scanning angle. First, ground-based mobile laser scanning (MLS) systems [65] should be taken into account for increasing the homogeneous distribution in the point cloud and thus enhancing a uniform field coverage. Second, unmanned aerial vehicles (UAVs), such as the recently introduced Riegl RiCOPTER [66], should be examined as a potential platform of a light-weight airborne laser scanning (ALS) systems. Promising results have already been achieved for measuring tree heights [67] or detecting pruning of individual stems [68] with UAV-based laser scanning. However, as examined in a comparative study for TLS and common plane-based ALS, the scanning angle and possible resolution influence the results [69] and thus have also to be taken into account for studies on UAV-based scanning systems.

In this study, TLS measurements were used to derive 3D information of points. As shown in other studies, captured intensity values could be used for qualitative analyses of the points, such as detecting single plants [70,71]. Whilst such analyses were not an object of this study they should be considered for further investigations. Moreover, full waveform analysis, commonly known from ALS, can simplify the distinction between laser returns on canopy and ground returns in TLS data [72,73]. The scanner used in this study however is not capable of capturing the full waveform.

The maps of plant height demonstrate the potential of the present approach for deriving plant height information on plot level in a very high resolution. The methodology of spatial plant height mapping can be scaled to field level, as long as the maximum range of the scanner is regarded and the point density is above the required minimum. As shown by Hämmerle and Höfle [63], the coverage of the field and attained mean heights are influenced by the point density. The approach of pixel-wise calculating plant height from TLS-derived CSMs has already shown good results at the field level for monitoring a maize field, about 80 m by 160 m in size [74] and a sugar beet field, about 300 m by 500 m in size [43] captured from four and eight scan positions, respectively. Further studies are necessary for determining crop- or case-specific minimum values for the point density. In this context, the used sensor and its maximum range influence the required number of scan positions.

Nevertheless, for this study, high coefficients of determination between averaged CSM-derived and manual measured plant heights validate the TLS measurements. For the absolute values, differences between the measurement methods have to be considered. Whereas for the manual measurement the heights of ten plants were averaged per plot, the CSM captured the entire crop surface. Consequently, differences in the mean heights occurred, which make precision analysis between TLS data and manual measurements infeasible. The precision of TLS measurements for agricultural applications is presumed from other studies [26,71]. It is important to note that a key advantage of the TLS data is that while plants for the manual measurements are subjectively selected, CSMs enable an objective assessment of spatially continuous plant height.

4.2. Biomass Estimation from Plant Height

Generally, PH performed well for the estimation of biomass in the pre-anthesis and the whole observed period. For dry biomass, PH was the best predictor for the whole observed period and similar good predictor as the best performing VIs for the pre-anthesis. However, PH performed far better for dry biomass than for fresh biomass, although these values are only distinguished by the water content of the sample. Thus, a possible explanation is the fact that the water content is not only influenced by the changing plant phenology across the growing season, but also by varying weather conditions. Moreover, during each day the available soil water and transpiration conditions vary. Hence, the amount of fresh biomass might vary more between the campaigns while the dry biomass is less influenced. Since PH is hardly affected by the water content of the plants, the varying water content in the fresh biomass adds noise to the BRM based on PH which results in lower R² values.

4.3. Biomass Estimation from Vegetation Indices

All VIs in this study have previously shown a relationship with biomass and LAI. Since the VIs use different bands within the spectral range, they were subdivided into three categories VIS VIs (RGBVI), VISNIR VIs (NDVI, RDVI, REIP), and NIR VIs (NRI, GnyLi). The VIs showed varying performances for the estimation of dry and fresh biomass, also depending on the regarded time frame of the growing season. Generally, the VIs within a category showed a similar behavior.

The saturation problem of the NDVI type VISNIR VIs was confirmed: Typically, crops reach 100% canopy cover around mid-vegetative phase. However, most crops continue to accumulate biomass and LAI afterwards. At a LAI of about 2.5–3, the absorbed amount of red light reaches a peak while the NIR scattering by leaves continues to increase. Thus, the ratio of NDVI type VISNIR VIs will only show slight changes [13]. In this study, the sensitivity thresholds were about 7.4 g/m² and 55 g/m² for dry and fresh biomass, respectively. Additionally, after heading the canopy de-greens due to flowering and fruit development (after BBCH 5, Table 2). This leads to an increased reflectance in the red part of the spectrum and thus, decreases values of the VISNIR VIs, while the biomass does not decrease. Herein, this discrepancy resulted in an inadequate model parameterization for the BRMs of the VISNIR VIs and poorer results for the whole observed period than for pre-anthesis.

A similar behavior was observable for the RGBVI. The inferior results might be explained by the fact that this VI does not take the reflectance in the NIR region into account, where most of the absorption features for biomass-related plant compounds are situated [75]. These results align well with the ones presented by Bendig et al. [54], where low correlations were found for the RGBVI with biomass after booting stage (BBCH 4, Table 2). In pre-anthesis, relationships of the RGBVI with dry and fresh biomass were similar. These results suggest that the RGBVI is mostly related with vegetation cover and not directly with biomass.

In contrast, NIRVIs, such as GnyLi and NRI, use bands only in the NIR and are thus not affected by the absorption in the red part of the spectrum, which could explain the overall more consistent and better performance of the NIR VIs, particularly after anthesis. A later saturation of these VIs aligns well with other studies [53,56]. Similarly, the REIP did not show any saturation effects in the pre-anthesis and yielded very good results for dry and fresh biomass. These findings can be explained by the major influence of the NIR bands that are not normalized as they are in the NDVI type VIs. Thus, the REIP saturated later than the VISNIR and VIS VIs. Nevertheless, across the whole observed period, the performance of the REIP also decreased due to saturation. The importance of the NIR domain for biomass estimation aligns with other studies [15,16,53,56] and should be further investigated. Similar to PH, the NIR VIs performed better for dry than for fresh biomass while the VISNIR VIs generally performed better with fresh biomass. This suggests that the VISNIR VIs respond more to the canopy water content and the related reflectance change in the NIR shoulder rather than directly to the biomass.

Overall, the results show that the NIR VIs perform best in the prediction of fresh and dry biomass. Moreover, the results indicate that the VIS and VISNIR VIs might not be directly related to biomass. However, no rigorous sensitivity analysis was carried out in this study but, as indicated by the results, such analyses should be carried out in the future.

In general, hyperspectral field measurements have been shown to be useful in earlier studies to estimate biomass [14,15,16,17]. However, VIs are prone to errors by illumination changes [76] and multiangular reflection effects [77]. So far, the influence of these effects on the estimation of plant parameters have not been comprehensively investigated and should be examined for evaluating the potential of VIs for plant parameter estimations. Moreover, ground-based spectrometer measurements are laborious and time-consuming. Automated platforms are under development in different fields of remote sensing to overcome this difficulty but they have not yet become standard. Kicherer et al. [78] developed a robotic platform for phenotyping grapevine based on automatic image acquisition. Results of a mobile multi-sensor phenotyping platform for phenotyping of winter wheat are presented by Kipp et al. [79]. Moreover, hyperspectral UAV-based systems showed promise [55,80,81,82,83]. Unfortunately, the promising NIR domain is currently not well covered by UAV sensing systems.

4.4. Biomass Estimation with Fused Models

Leaves make up a major part of the biomass, and VIs related to biomass are often also responsive to LAI [13,84]. Thus, it was assumed that the spectral information would complement the PH information by adding information about the canopy density and cover.

As described above, PH and VIs showed varying performances in the estimation of fresh and dry biomass and for pre-anthesis or the whole observed period. For dry biomass in pre-anthesis, the NIR VIs performed slightly better than PH. Here, the fusion with all VIs improved the predictability, whereby the NIR and VISNIR VIs yielded the best results. This can be explained by the sensitivity of the VIs to the vegetation cover in early growth stages. For the whole observed period, PH clearly outperformed the VIs in the multivariate BRMs and only the fusion with the NIR VIs increased the predictability slightly compared to PH alone. For the VIS and VISNIR VIs, the above described saturation effects might have counteracted the positive effect of the vegetation cover estimation in the early growth stages. Additionally, for pre-anthesis and across the whole observed period, the multivariate BRMs performed similarly regardless which VI was used. This indicates that most of the prediction power can be accounted to PH.

For fresh biomass across the whole observed period, the NIR VIs performed best, followed by PH. Although the VISNIR VIs did not perform well in the bivariate BRMs, they could improve the results when fused with PH. As described above, VISNIR VIs respond to the water content. Thus, they might have complement the PH information for the estimation of fresh biomass. Still, only a slight improvement was achieved with the fused models compared to the NIR VIs alone and overall, the results of multivariate BRMs with different VIs differed only slightly.

In pre-anthesis, only the NDVI and RGBVI performed poorer than PH while the REIP performed best for the fresh biomass. In combination with PH, the results of the NDVI and RDVI were improved the most, while the latter one also achieved the best results of all fused models. For the NIR VIs and REIP none or only very minor improvements were achieved and as for the whole observed period, the water was important because it influences the reflectance in the NIR. Additionally, the VIs correspond to vegetation cover in the early growth stages. Thus, in pre-anthesis already the VIs performed well and PH only rarely contributed to the prediction power. Only the RGBVI, NDVI, and RDVI might have carried complementary information to the PH.

In this study, the NIR VIs showed the overall best performance of the VIs and seemed to carry similar information as PH. Overall, PH and NIR VIs showed the best potential for biomass estimation as individual and fused estimators. This aligns with a recent study by Marshall and Thenkabail [15], in which they have shown the importance of PH and the NIR domain for fresh biomass estimations. The VISNIR VIs seemed to be influenced by the water content and their performance strongly depended on the regarded time frame of the growing season. Although, no comprehensive sensitivity analysis was carried out, these findings align well with other studies [56,85]. Further studies are needed to investigate the influence of the growing stage on the estimation, and whether estimators, which have been found as suitable in across growth stage estimations, are suitable for estimation at individual growth stages. Such in-season estimations are particularly important for applications in precision agriculture. Additionally, in this study VIs known for estimating biomass from hyperspectral data were used. Thus, the full potential of the fusion of 3D spatial and spectral data may not have been explored. Future studies should investigate whether other parts of the spectral range complement PH information better.

Overall, this study demonstrated the strength of bivariate BRMs based on PH and NIR VIs for estimating biomass, with only slight improvements achievable through multivariate models. In contrast, the weak performances of the VIS and VISNIR VIs as individual estimators were compensated through the fusion with PH. However, statements have to be limited, since the models indicated that PH contributed the most to the prediction power. In this context, it has to be noted that neither linear nor exponential models reflected the relation between estimators and biomass perfectly and thus more complex functions have to be considered, which might take the benefits of VIs, like sensitivity to water content, better in to account.

For practical applications the benefit of the fused models might be outweighed by the expenses to deploy two different systems. Referring to this, limitations through the attainable spatial and temporal resolution of each system have to be regarded. As already mentioned, TLS measurements can be scaled up to larger fields, as long as a sufficiently point density can be achieved, which has to be determined crop- and case-specific in further studies. Apart from that, laser scanning appeared as powerful tool for the non-destructive and objective assessment of spatially resolved plant height data. Statements about the accuracy of the measured plant heights are hardly possible due to the already mentioned different spatial resolution of the plant height measurements, however the averaged difference of 0.05 m between TLS-derived and manual measured plant heights corroborate the results (Table A1). A main benefit of the field spectrometer measurements is the high credibility of the acquired spectral data, based on a large number of former studies, however the dependence on solar radiation and the small numbers of measurements per regarded spatial area, herein per plot, are the main disadvantages. Consequently, systems are required which are capable to assess larger areas in less time with the same accuracy of the results. Ideally, spatial and spectral information should be acquired directly through one sensor. For example, recently developed sensing systems and techniques allow to create hyperspectral point clouds [86] and hyperspectral digital surface models [83] with only one sensor and thus, derive 3D spatial and hyperspectral information at the same time. Thus, it can be expected that 3D hyperspectral information will become increasingly available and combined analysis approaches should be further developed.

5. Conclusions and Outlook

Continuously conducting a field experiment with different barley cultivars and the related TLS, field spectrometer, and manual measurements enabled the acquisition of an extensive data set. High R² values up to 0.89, between TLS-derived and manual measured plant heights verified the applicability of the presented approach for a pixel-wise calculation of plant height (PH) from high resolution crop surface models (CSMs). Six established vegetation indices (VIs) were used to extract information from the hyperspectral data. Based on PH and VIs, bivariate and multivariate biomass regression models (BRMs) were established, with varying performances. Whereas PH was supported as strong estimator in the bivariate models (R² up to 0.85), VIs showed highly different results (R²: 0.07–0.87). The multivariate models yielded improvements in some cases (R² up to 0.89), however in most cases PH had the greatest contribution to the prediction power.

Different models appeared best suitable for dry or fresh biomass estimations, also depending on the regarded time frame of the growing season, but in all cases exponential models performed better than the linear ones: For dry biomass, the bivariate BRM with PH showed the best results for the whole observed period (R² = 0.85), whereas for the pre-anthesis the REIP and the near-infrared (NIR) VIs GnyLi and NRI showed slightly better results than PH (R² = 0.86, 0.87). Multivariate BRMs from PH and one VI slightly improved the R² values compared to the bivariate BRMs in some cases. For fresh biomass, the bivariate BRMs of the NIR VIs showed the best results for the whole observed period (both R² = 0.77). For pre-anthesis, the REIP (R² = 0.82) showed slightly better results that the NIR VIs (both R² = 0.79). The multivariate BRM could slightly improve the results in some cases. Additionally, it can be noted that also weakly performing VIs, such as the NDVI or RGBVI, improved the estimations slightly when fused with PH in the multivariate BRMs, both for fresh and dry biomass. These results suggest that specific models should be chosen for specific applications, and a fusion of PH and VIs does not always substantially improve the results. Additionally, when PH and VIs are fused, the choice of the VI does not seem critical in all cases.

Altogether, it should be noted that the presented results are a first step towards the fusion of remotely sensed 3D spatial and spectral data for a precise and non-destructive estimation of crop biomass. Other ways of data fusion may further increase the prediction power. Further studies are also necessary to investigate differences between the years, cultivars, and fertilizer treatments. Moreover, as already mentioned, in-season biomass estimations are important for precision agriculture. Therefore models should be established based on data sets from only one campaign to investigate the potential for timely monitoring and in-season estimations. Accurate and rapidly ascertainable estimations in a high spatial resolution during the growing season could support spatially resolved nitrogen nutrition index calculations. Thereby in-field variations can be considered for optimizing fertilizer application and shrinking the gap between potential and current yield. The fusion of 3D spatial and spectral data might improve such calculations as weaknesses and limitations of one estimator might be compensated through the other one.

With regard to the application in the field, the usability of new platforms should be further investigated. UAV-based light-weight ALS systems reveal potential for vegetation mapping. Futhermore, new technologies like hyperspectral snapshot camera systems which enable the derivation of 3D spatial and hyperspectral information at the same time carry great potential for agricultural applications. Combined with estimation models based on structural and spectral and information, such approaches can become a powerful tool for applications in precision agriculture and biomass monitoring.