Predicting Forage Quality of Grasslands Using UAV-Borne Imaging Spectroscopy

: The timely knowledge of forage quality of grasslands is vital for matching the demands in animal feeding. Remote sensing (RS) is a promising tool for estimating ﬁeld-scale forage quality compared with traditional methods, which usually do not provide equally detailed information. However, the applicability of RS prediction models depends on the variability of the underlying calibration data, which can be brought about by the inclusion of a multitude of grassland types and management practices in the model development. Major aims of this study were (i) to build forage quality estimation models for multiple grassland types based on an unmanned aerial vehicle (UAV)-borne imaging spectroscopy and (ii) to generate forage quality distribution maps using the best models obtained. The study examined data from eight grasslands in northern Hesse, Germany, which largely di ﬀ ered in terms of vegetation type and cutting regime. The UAV with a hyperspectral camera on board was utilised to acquire spectral images from the grasslands, and crude protein (CP) and acid detergent ﬁbre (ADF) concentration of the forage was assessed at each cut. Five predictive modelling regression algorithms were applied to develop quality estimation models. Further, grassland forage quality distribution maps were created using the best models developed. The normalised spectral reﬂectance data showed the strongest relationship with both CP and ADF concentration. From all predictive algorithms, support vector regression provided the highest precision and accuracy for CP estimation (median normalised root mean square error prediction (nRMSE p ) = 10.6%), while cubist regression model proved best for ADF estimation (median nRMSE p = 13.4%). The maps generated for both CP and ADF showed a distinct spatial variation in forage quality values for the di ﬀ erent grasslands and cutting regimes. Overall, the results disclose that UAV-borne imaging spectroscopy, in combination with predictive modelling, provides a promising tool for accurate forage quality estimation of multiple grasslands.


Introduction
In Europe, approximately 30% to 35% of the agricultural area consists of grasslands [1]. Mainly permanent grasslands are incredibly variable in species composition, biodiversity, and management practices, as well as in productivity [2]. Food provisions as forage for ruminants and herbivores and as biomass substrate for energy production are the most comprehensive ecosystem services from grasslands. There exists a multitude of destructive and non-destructive methods to measure or estimate the production and quality of the forage [3][4][5]. Timely and accurate information about forage quality is essential to match available feed on offer with animals' demand.
Gaussian processing regression (GPR), support vector regression (SVR), and cubist regression (CBR), that have not been examined with spectral data to estimate forage quality.
Typically, the models developed to estimate forage quality using spectral data are restricted to single grassland types. Consequently, the transferability of the models to other grassland types is limited due to the low variability of the underlying training data. Therefore, a general model that can estimate forage quality parameters, irrespective of the grassland type, would be a preferable solution. Thus, we hypothesise that the usage of spectral data over different grasslands with predictive modelling algorithms can build models to estimate major forage quality parameters regardless of the grassland type or management practice. The overall objective of this study is to develop comprehensive predictive models to estimate CP and ADF concentrations (from now on referred to as CP and ADF) on the field level, using UAV-borne imaging spectroscopy data from different grasslands. Specifically, to build models that can predict values regardless of the grassland type and to examine the effect of grassland type for the model performances. With this in mind, the specific sub-objectives are as follows: • To understand the relationship between forage quality parameters (CP and ADF) and spectral reflectance; • To develop models for CP and ADF estimation from imaging spectroscopy data and to evaluate models; • To create field-level CP and ADF maps for grasslands and describe the spatial and temporal variation of forage quality.

Study Area
The study was carried out on eight grasslands with different management practice and species composition. Table 1 gives details of the grasslands, and Figure 1 shows their locations on the map. Six grasslands were in nature protection areas, where no fertilisation was applied, and they were mowed only once per year. Two of them were mountain hay meadows (MHM), and another two were Nardus stricta grasslands (NSG). MHML and NSGL were MHM and NSG grasslands, which were substantially invaded by the neophyte Lupinus polyphyllus. The lowland hay meadow (LHM) is an extensively utilised grassland located on the Werra riverbank, and it was mowed two times per year. The intensively managed grassland (IMG) was fertilised and was harvested three times per year. All grasslands were in northern Hesse, Germany. MHM1, NSG1, LHM, and IMG grasslands were in Werra-Meißner district in Hesse (EPSG 4326: 9.9 • N, 51. A 25 m by 50 m rectangle plot (1250 m 2 ) was selected as the study plot in grasslands located in Werra-Meißner district. In each study plot, twenty 1 m 2 subplots were randomly distributed and selected for destructive biomass samplings. The random sub-plot generation was done in QGIS, with each subplot having a minimum distance of 5 m to the nearest subplot.
The plots in the biosphere reserve Rhön were not allowed to be mown before 15 July and, to replicate different vegetation maturity stages, different sampling design was employed. From the grasslands in biosphere reserve Rhön, a 50 m by 30 m rectangle plot (1500 m 2 ) was chosen as the study plot. A total of fifteen small plots (8 m by 8 m) were established in each study plot, and they were in a regular grid (five rows and three columns). Three 1 m 2 sampling plots were picked from each small plot as the subplots for destructive biomass sampling. A 25 m by 50 m rectangle plot (1250 m 2 ) was selected as the study plot in grasslands located in Werra-Meißner district. In each study plot, twenty 1 m 2 subplots were randomly distributed and selected for destructive biomass samplings. The random sub-plot generation was done in QGIS, with each subplot having a minimum distance of 5 m to the nearest subplot.
The plots in the biosphere reserve Rhön were not allowed to be mown before 15 July and, to replicate different vegetation maturity stages, different sampling design was employed. From the  The Cubert Hyperspectral Firefleye S185 SE (Cubert GmbH, (Ulm, Germany) www.cubert-gmbh. de) snapshot camera (12 mm focal length) is an imaging spectroscopy sensor that acquires spectral images between 450 and 998 nm. It is a 2D imager with a multi-point spectrometer. The spectral resolution of the sensor (full width at half maximum) is 4.8 nm at 450 nm and 25.6 nm at 850 nm. The camera records a total of 138 spectral bands with a 4 nm sampling interval, but the manufacturer recommended to use only 126 bands between 450 and 950 nm. Each spectral band image is 50 by 50 pixel, and radiometric resolution is 12 bit (0-4096 digital numbers). In addition to the spectral bands, a grayscale panchromatic image is also recorded with 1000 by 990 pixels in size [36,37]. The spatial resolution of the spectral image is 20 cm at 20 m flying height. The captured image from the camera is an image cube that has 126 bands.
The Cubert camera was attached to the UAV (RTK X8 Hyperspectral Mapping, Cubert GmbH, www.cubert-gmbh.de). It is a co-axial multi-rotor UAV equipped with real-time kinematic (RTK) global navigation satellite system (GNSS). It can fly automatically guided with less than 1 m position accuracy. The maximum payload is 2.5 kg, and maximum flight weight is 10 kg. The camera is affixed to a gimble on the UAV, which compensates for pitch and roll movement during the flight, and keeps the camera always nadir looking [38]. Figure 2 shows the UAV and the attached Cubert camera. The Cubert Hyperspectral Firefleye S185 SE (Cubert GmbH, (Ulm, Germany) www.cubertgmbh.de) snapshot camera (12 mm focal length) is an imaging spectroscopy sensor that acquires spectral images between 450 and 998 nm. It is a 2D imager with a multi-point spectrometer. The spectral resolution of the sensor (full width at half maximum) is 4.8 nm at 450 nm and 25.6 nm at 850 nm. The camera records a total of 138 spectral bands with a 4 nm sampling interval, but the manufacturer recommended to use only 126 bands between 450 and 950 nm. Each spectral band image is 50 by 50 pixel, and radiometric resolution is 12 bit (0-4096 digital numbers). In addition to the spectral bands, a grayscale panchromatic image is also recorded with 1000 by 990 pixels in size [36,37]. The spatial resolution of the spectral image is 20 cm at 20 m flying height. The captured image from the camera is an image cube that has 126 bands.
The Cubert camera was attached to the UAV (RTK X8 Hyperspectral Mapping, Cubert GmbH, www.cubert-gmbh.de). It is a co-axial multi-rotor UAV equipped with real-time kinematic (RTK) global navigation satellite system (GNSS). It can fly automatically guided with less than 1 m position accuracy. The maximum payload is 2.5 kg, and maximum flight weight is 10 kg. The camera is affixed to a gimble on the UAV, which compensates for pitch and roll movement during the flight, and keeps the camera always nadir looking [38]. Figure 2 shows the UAV and the attached Cubert camera.

Spectral Images and Forage Quality Data
UAV-borne imaging spectroscopy data collection and destructive biomass sampling were carried out on every grassland immediately prior to each cut between May and September in 2018 ( Table 2). Table 2. Details on samplings of the grasslands investigated in the study.

Spectral Images and Forage Quality Data
UAV-borne imaging spectroscopy data collection and destructive biomass sampling were carried out on every grassland immediately prior to each cut between May and September in 2018 (Table 2). Third cut 1 August 20 † † 1 m 2 plots. ‡ at each cut, three 1 m 2 subplots inside the 64 m 2 plots were sampled, which gives a total of fifteen samples (3 × 5) at each cut for one grassland. However, in the lab analysis, samples from the three sub-plots were mixed and concentration data were determined (five quality samples per cut in one grassland).
Spectral images were acquired between 10:00 and 14:00 under clear-sky conditions. Before the UAV flight, four corners of the study area were staked out using the Leica RTK GNSS (Leica Geosystem, www.leica-geosystems.com) and six black-white 1 m 2 ground control points (GCP) were distributed around the study plot. The coordinates of the GCPs were also measured using the RTK GNSS. The UAV flight mission was done with guided auto-pilot system in parallel flight lines at a flying height of 20 m and a horizontal flying speed of 5 ms −1 . Image acquisition points were calculated based on 80% image overlap in both directions. In each image capturing point, the UAV hovered for 1 s to minimise the shaking effect. The UAV automatically triggered the Cubert camera on each correct position.
Before the flight, the radiometric calibration for the Cubert camera was done using a spectral reference image on the white calibration panel (95% reflectance Zenith Lite) and a dark current image by closing the lens using a lens cap. The relative spectral reflectance value for each image pixel was calculated according to Equation (1) during image acquisition.
where i, j, and λ represent image row number, column number, and wavelength, respectively; R is the relative spectral reflectance; and DC, WC, and DN are the digital numbers of the dark current image, the white calibration panel image, and the raw image, respectively [39]. After spectral data collection, grass biomass was clipped on the identified 1 m 2 subplots at a stubble height of 5 cm. The fresh biomass was weighted in the field, and total biomass in each plot was divided into two separate samples for dry biomass and forage quality analysis. The samples for the quality analysis were dried at 65 • C for 48 hours. Afterwards, the dried samples were ground for 1 mm uniform particles with a Foss CT 193 Cyclotec mill (FOSS, www.fossanalytics.com). Subsequently, the dry matter and the ash content of the ground samples were determined, and the N and ADF were determined.
The ADF was determined using the ANKOM 200 Fibre Analyser (ANKOM Technology, www. ankom.com). The samples were cooked in the analyser at 100 • C for one hour in strong acid detergent (Cetyltrimethylammonium bromid in 1 N H 2 SO 4 ). Afterwards, the cooked samples were washed, dried, and weighed to determine the ADF. The N concentration was measured using a Vario Max CHN elemental analyser (Elementar Analysensysteme GmbH, (Langenselbold, Germany) www.elementar.de), which employs the Dumas method in analytical chemistry. CP was computed multiplying the N concentration by 6.25. CP and ADF were determined for a total of one hundred and ninety-four (194) samples in all grasslands (Table 2).

From Image Cubes to Digital Ortho-Mosaic
The following procedure was applied in each image acquisition to generate a mosaic from Cubert image cubes ( Figure 3). First, the image cubes were exported as TIFF files using the CubeExport DOS (Cubert GmbH, www.cubert-gmbh.de) programme. It produced a stacked spectral image (50 × 50 pixels) and a grayscale image (1000 × 1000 pixels). Afterwards, the spectral images were disaggregated (nearest neighbour resampling) 20 times and stacked with the grayscale image using raster package in R programming language [40,41]. Then, the pre-processed stacked images were input to the Agisoft PhotoScan Professional version 1.4.1 (64 bit) software (Agisoft LLC, (St. Petersburg, Russia) www.agisoft.com) [42]. The PhotoScan software was used to create the ortho-mosaic from single images. The measured GCP coordinates were fed to the PhotoScan during the mosaic generation. The geo-referenced mosaic was then corrected for elevation to produce an ortho-mosaic with 127 layers (grayscale band + 126 spectral bands). The mean spectral reflectance data for each 1 m 2 subplot were extracted from the ortho-mosaics of each sampling date. Due to noise in the spectral bands between 450 and 478 nm (8 bands), only 118 spectral bands (482-950 nm) were utilised in this study.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 23 The following procedure was applied in each image acquisition to generate a mosaic from Cubert image cubes ( Figure 3). First, the image cubes were exported as TIFF files using the CubeExport DOS (Cubert GmbH, www.cubert-gmbh.de) programme. It produced a stacked spectral image (50 × 50 pixels) and a grayscale image (1000 × 1000 pixels). Afterwards, the spectral images were disaggregated (nearest neighbour resampling) 20 times and stacked with the grayscale image using raster package in R programming language [40,41]. Then, the pre-processed stacked images were input to the Agisoft PhotoScan Professional version 1.4.1 (64 bit) software (Agisoft LLC, (St. Petersburg, Russia) www.agisoft.com) [42]. The PhotoScan software was used to create the orthomosaic from single images. The measured GCP coordinates were fed to the PhotoScan during the mosaic generation. The geo-referenced mosaic was then corrected for elevation to produce an orthomosaic with 127 layers (grayscale band + 126 spectral bands). The mean spectral reflectance data for each 1 m 2 subplot were extracted from the ortho-mosaics of each sampling date. Due to noise in the spectral bands between 450 and 478 nm (8 bands), only 118 spectral bands (482-950 nm) were utilised in this study.

The Relationship between Reflectance Data and Forage Quality
Correlation analysis was carried out to understand the relationship between CP and ADF in the grassland biomass and the reflectance data. First, spectral reflectance data were normalised using vector normalisation based on Equation (2) [43], where xi is the spectral vector for i = 1, 2,…, n. The normalised reflectance data was transformed into the first derivative reflectance [44] and continuum removal band depth [45].
Pearson correlation coefficients (r) between various transformed reflectance (actual, normalised, derivative, and continuum removal) data and forage quality data were computed. The optimal transformation method was identified based on the highest average correlation coefficient value.

The Relationship between Reflectance Data and Forage Quality
Correlation analysis was carried out to understand the relationship between CP and ADF in the grassland biomass and the reflectance data. First, spectral reflectance data were normalised using vector normalisation based on Equation (2) [43], where x i is the spectral vector for i = 1, 2, . . . , n. The normalised reflectance data was transformed into the first derivative reflectance [44] and continuum removal band depth [45].
Pearson correlation coefficients (r) between various transformed reflectance (actual, normalised, derivative, and continuum removal) data and forage quality data were computed. The optimal transformation method was identified based on the highest average correlation coefficient value. Simple linear regression models were created between each band from the identified transformed reflectance data set and forage quality data.
Further, all possible variants of SR and NDSI (Equations (3) and (4)) were derived from the reflectance data to check for the relationship with CP and ADF by simple linear regression modelling. The resulting models were evaluated based on the adjusted coefficient of determination (adj.R 2 ) value.
where λ 1 , λ 2 are two wavelengths and ρ is the reflectance at the wavelength.

Forage Quality Modelling with Full Spectral Data
In conformity with the results from the simple linear regression analysis, neither single wavelengths nor SR or NDSI provided strong relationships with CP or ADF (data are shown in Results in Section 3.3) Consequently, full spectral reflectance data were used instead of selected spectral features. Five predictive modelling algorithms were utilised, namely, partial least squares regression (PLSR), random forest regression (RFR), Gaussian processes regression (GPR), support vector regression (SVR), and cubist regression (CBR). A short description of each regression algorithm with references is given in Table 3. Table 3. Summary of the predictive modelling algorithms used in the study.

Algorithm
Description Reference PLSR Partial least squares regression builds a linear regression model on the data projected in a space, based on nonlinear iterative partial least squares [46] GPR Gaussian processes regression finds a regression solution based on a probabilistic approach [47,48] RFR Random forest is an ensemble method consisting of decision trees and bagging [49] SVR Support vector regression builds linear regression models in a high-dimensional feature space [50][51][52] CBR Cubist is a rule-based regression technique with boosting functionality [53][54][55] Predictive modelling can be either supervised or unsupervised learning methods [56]. Thus, a supervised learning approach was employed in this study. First, the data set was split into two portions, representing the training and test dataset. Then, the model was trained using only the train data, while the test data were held out. After the model was trained, the model performance was evaluated using the test data (hold-out method) [57]. However, one random split of the data in the hold-out method results in only one model and the model performance may differ significantly among different split sets from the data [56]. Developing separate models for different data splits could explain the variability in model performance within the given data subsets. Therefore, a simulation was done 100 times and each replicate consisted of separate training (80% data) and test (20% data) data sets. Data set split was based on the stratified random sampling method with grassland field as strata. Then, 100 separate models using 100 training data sets were built for each regression algorithm using the R package "classification and regression training (caret)" [40,56,58]. In the model training phase, the repeated cross-validation (five folds; three repeats) technique was applied to tune the models' hyper-parameters (ncomp in PLSR; sigma in GPR; mtry in RFR; sigma & cost in SVR; committees & neighbours in CBR).
Hyper-parameters are model configuration parameters that optimise the performance of the predictive modelling algorithm [57]. Finally, model testing was done using the corresponding 100 different test data. To evaluate the prediction performances of the model, the squared correlation coefficient between observed and prediction R 2 p [59], root mean square error prediction (RMSE p ) (Equation (5)), and normalised RMSE p based on observation range (nRMSE p ) (Equation (6)) [60] were calculated. The overall workflow is shown in Figure 4.
where y is the observed value,ŷ is the predicted value, y max is the maximum observed value, y min is the minimum observed value, and n is the number of observations. nRMSE p = RMSE p y max -y min ×100, where y is the observed value, y is the predicted value, ymax is the maximum observed value, ymin is the minimum observed value, and n is the number of observations. On the basis of the optimal regression algorithms and hyper-parameter values, two models were calibrated using the 100% data with the repeated cross-validation technique (five folds; three repeats) as the final predictive models for CP and ADF. The performance of final models was evaluated using R 2 cv [59], RMSEcv (Equation (5)), and nRMSEcv (Equation (6)).

Forage Quality Prediction Maps
Model predictions were utilised to create forage quality distribution maps. For this purpose, the ortho-mosaic was averaged to 1 m 2 pixel to match the ground sampling distance of the models, and for every pixel, CP and ADF were predicted for each harvest in each grassland. The predicted forage quality maps were compared visually for different cuts in the same grassland and examined for spatial and temporal variations. On the basis of the optimal regression algorithms and hyper-parameter values, two models were calibrated using the 100% data with the repeated cross-validation technique (five folds; three repeats) as the final predictive models for CP and ADF. The performance of final models was evaluated using R 2 cv [59], RMSE cv (Equation (5)), and nRMSE cv (Equation (6)).

Forage Quality Prediction Maps
Model predictions were utilised to create forage quality distribution maps. For this purpose, the ortho-mosaic was averaged to 1 m 2 pixel to match the ground sampling distance of the models, and for every pixel, CP and ADF were predicted for each harvest in each grassland. The predicted forage quality maps were compared visually for different cuts in the same grassland and examined for spatial and temporal variations.

Forage Quality Data
Crude protein concentration varied between 5.1 %DM and 23.3 %DM, while ADF varied between 22.5 %DM and 38.5 %DM (Table 4). Forage from the intensively managed grassland (IGM) had the highest average CP and the lowest average ADF ( Figure 5). However, forage from the two grasslands invaded by Lupinus polyphyllus (MHML, NSGL) contained higher CP than non-invaded grasslands. Further, CP and ADF were inversely correlated (correlation coefficient −0.78).

Imaging Spectroscopy Data.
Normalised mean reflectance values in the visible region (500-700 nm) obtained for every 1 m 2 sampling plot were lower for samples with higher CP values, along with higher values in the nearinfrared region (750-900 nm) ( Figure 6). A similar pattern was found for ADF data. However, there

Imaging Spectroscopy Data
Normalised mean reflectance values in the visible region (500-700 nm) obtained for every 1 m 2 sampling plot were lower for samples with higher CP values, along with higher values in the near-infrared region (750-900 nm) ( Figure 6). A similar pattern was found for ADF data. However, there were some reflectance curves that did not follow the natural vegetation spectral curve, as they probably contained mixed vegetation signals.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 23 were some reflectance curves that did not follow the natural vegetation spectral curve, as they probably contained mixed vegetation signals.

Linear Modelling with Individual Bands and Spectral Indices (Spectral Features)
Correlograms for the different reflectance data sets (original, normalised, derivatives, and continuum removal band depth) are shown in Figure 7. On the basis of correlation values, normalised reflectance data obtained the highest average r-value from all wavelengths (0.26 for both CP and ADF) and highest cumulative sum of r-value (30.5 for both CP and ADF). Therefore, the normalised reflectance data were used for further statistical analysis in this study. For CP, the normalised reflectance value obtained negative r-values in the visible region, but these became positive in the near-infrared region. The correlogram of normalised reflectance for ADF showed an opposite pattern compared with that of CP.

Linear Modelling with Individual Bands and Spectral Indices (Spectral Features)
Correlograms for the different reflectance data sets (original, normalised, derivatives, and continuum removal band depth) are shown in Figure 7. On the basis of correlation values, normalised reflectance data obtained the highest average r-value from all wavelengths (0.26 for both CP and ADF) and highest cumulative sum of r-value (30.5 for both CP and ADF). Therefore, the normalised reflectance data were used for further statistical analysis in this study. For CP, the normalised reflectance value obtained negative r-values in the visible region, but these became positive in the near-infrared region. The correlogram of normalised reflectance for ADF showed an opposite pattern compared with that of CP.
With simple linear regression modelling, a weak relationship was found for both CP and ADF against individual spectral bands ( Table 5). The best performing model for CP was found with 718 nm (adj.R 2 = 0.33) and for ADF with 794 nm (adj.R 2 = 0.23). Overall, the model performance was better for CP in the visible region and, vice-versa, for ADF in the near-infrared region.
reflectance data obtained the highest average r-value from all wavelengths (0.26 for both CP and ADF) and highest cumulative sum of r-value (30.5 for both CP and ADF). Therefore, the normalised reflectance data were used for further statistical analysis in this study. For CP, the normalised reflectance value obtained negative r-values in the visible region, but these became positive in the near-infrared region. The correlogram of normalised reflectance for ADF showed an opposite pattern compared with that of CP.   Linear regression models with NDSI and SR performed better than single-band models (Table 5). Models with the two indices performed similarly, and the best results were found for similar band combinations for both forage quality values. For CP estimation, the best NDSI and SR models obtained adj.R 2 values of 0.42 and 0.40, respectively. Similarly, for ADF, the best NDSI and SR models obtained adj.R 2 values of 0.34 and 0.33, respectively.

The Best Predictive Modelling Algorithm
On the basis of 100 different model runs with random train and test data sets, the model with the highest accuracy (lowest median RMSE p and lowest median nRMSE p ), the highest precision (highest R 2 p ), and the highest stability (lowest standard deviation of RMSE p ) was identified as the best-performing model ( Table 6 and Figure 8). Accordingly, the SVR model (median RMSE p = 1.9 %DM; median nRMSE p = 10.6%; median R 2 p = 0.79; SD RMSE p = 0.29 %DM) was the best model for CP estimation, whereas for ADF, the CBR model (median RMSE p = 2.2 %DM; median nRMSE p = 13.4%; median R 2 p = 0.56; SD RMSE p = 0.23 %DM) was the optimal model. PLSR was the least performing model type among all the predictive algorithm models. For CP estimation, nRMSE p of the SVR model varied between 7.0% and 14.5% and nRMSE p varied from 6.5% to 16.4% of the CBR model. Comparably, nRMSE p of ADF models varied from 10.7% to 18.3% and from 10.5% to 16.7% for the SVR and CBR Remote Sens. 2020, 12, 126 14 of 23 models, respectively. Further, the precision of (median R 2 p ) CP models was larger than 0.73, except for the PLSR model. However, the precision of (median R 2 p ) ADF models was lower than for CP models.  The plots of fit for the best-performing models show the model fit across all grasslands ( Figure  9). Overall, prediction accuracy tended to be lower at higher levels of CP, whereas for ADF, accuracy was consistent across the whole range of values observed. The plots of fit for the best-performing models show the model fit across all grasslands (Figure 9). Overall, prediction accuracy tended to be lower at higher levels of CP, whereas for ADF, accuracy was consistent across the whole range of values observed.  Apart from the determination of the optimal predictive modelling algorithm, the best hyper-parameter values for each algorithm were identified. In the SVR model for CP estimation, the best sigma value was 0.035, and the optimum cost parameter was 20. In the CBR model for ADF, the optimal committees and neighbours values were 91 and 9, respectively (Table 7).

Final Models
With SVR and the CBR predictive algorithms being identified as the best algorithms to estimate CP and ADF, repeated cross-validations were performed using the complete data set. The SVR model resulted in a nRMSE cv of 9.6% with R 2 cv = 0.81 for CP estimation, while for ADF estimation, the CBR model had a nRMSE cv of 13.0% and a R 2 cv of 0.60. The errors of the final models were found between the errors from 100 different models in the model training and testing phase.

Forage Quality Maps
The final models (SVR for CP and CBR for ADF) were applied to the mosaic images to generate forage quality prediction maps exemplary for the IMG and MHML grasslands (Figures 10 and 11). The produced maps show observed spatial and temporal variability for both quality parameters. The maps of the IMG grassland ( Figure 10) indicate an increase in CP from the first to the second cut and a decrease from second to the third cut, which agrees with the values of the reference data ( Figure 5). The maps for ADF (Figure 10g-i) indicate an opposite pattern compared with CP (Figure 10d-f). Both CP and ADF maps for the second cut (Figure 10e,h, respectively) indicate a considerably higher spatial variation compared with the first and third cut. The ADF map of the third cut (Figure 10i) shows certain areas (towards the west-southwest corner) with somewhat lower values (purple colour patches) compared with the rest of the field.   The extracted CP and ADF values from each cutting date of the MHML grassland were plotted against each other to understand the temporal variation of the relationship between the two forage quality parameters (Figure 12). The strength of the negative linear relationship between the two parameters decreased when the cutting date was delayed. Further, the distribution of CP values decreased with increasing duration of growth, but the distribution of ADF values did not change significantly. The true colour images of the MHML grassland show a unique chessboard pattern at the later cutting dates (Figure 11). While on the 13th June, the whole area had the same maturity stage (the first cut was taken at 15 June on a third of all plots), the five brownish areas in the picture from 27 June represent the plots that were regrown since the first cutting date. Brown colour shades further increase in the picture from 11 July. The forage quality maps for the last two dates were masked using small plot boundaries (64 m 2 ). The CP map from 13 June shows higher CP value patches (yellow patches) at the northern region of the plot and lower ADF values (purple spots) for corresponding areas are shown in the ADF map from the same date. The most matured grass plots were available on 11 July, where the lowest CP values are shown (dark purple plots). Those patterns are less pronounced in the ADF map of 11 July. Colours on the map from 27 June indicate that the biomass on recently cut areas has more significant ADF concentrations than that from areas with continued growth, and thus an advanced plant maturity.
The extracted CP and ADF values from each cutting date of the MHML grassland were plotted against each other to understand the temporal variation of the relationship between the two forage quality parameters ( Figure 12). The strength of the negative linear relationship between the two parameters decreased when the cutting date was delayed. Further, the distribution of CP values decreased with increasing duration of growth, but the distribution of ADF values did not change significantly.

Discussion
Correlograms were used to understand the relationship between grassland CP and ADF and spectral reflectance (Figure 7). For both quality parameters, normalised spectra indicated a stronger relationship with forage quality parameters than the original spectra or other spectra transformations (derivatives and continuum removal band depth). The normalisation process reduced the temporal noise of the data, which strengthened the relationship between spectral features and forage quality. Further, two correlograms from CP and ADF showed an inverse pattern, confirmed by the negative correlation between them. Nevertheless, both forage quality parameters had a stronger relationship in the red-edge and the near-infrared regions of the spectrum.
Linear regression models with spectral features (single waveband, NDSI, SR) obtained a maximum adj.R 2 of 0.42 for CP and 0.34 for ADF. The best correlated NDSIs comprised wavebands in the chlorophyll absorption features in the visible domain of the spectrum (470-518 nm and 550-750 nm), which relates to N and other biochemicals in green canopies [14]. These results paralleled with the results from previous studies that come from the same geographical region [15,18]. The study from Safari et al. [32] obtained linear regression models with NDSI with R 2 of 0.58 and 0.49 for CP and ADF, respectively. Moreover, the results from Biewer et al. [15] reported R 2 values of 0.33 and 0.13 R 2 for CP and ADF models with SR, respectively. Nonetheless, the results from the studies mentioned above were obtained with field spectroscopy data, which had a more comprehensive spectral range (305-1700 nm) than the camera utilised in the present study.
Several predictive modelling algorithms were tested to identify the best algorithms to estimate CP and ADF from the full spectral data, as no single consistent algorithm was shown to surpass all the given circumstances every time in a study by the authors of [61]. Moreover, model consistency was evaluated by training and testing 100 models using 100 different random train and test data sets, which allowed to disclose the model performance irrespective of the calibration data set. Except for the PLSR, other tested predictive algorithms (GPR, RFR, SVR, and CBR) proved promising for the underlying data. The SVR and CBR models for CP and ADF estimation showed the maximum precision (highest R 2 p) and prediction accuracy (lowest nRMSEp) followed by the RFR, GPR, and PLSR models, respectively. SVR and CBR predictive modelling algorithms were previously utilised

Discussion
Correlograms were used to understand the relationship between grassland CP and ADF and spectral reflectance (Figure 7). For both quality parameters, normalised spectra indicated a stronger relationship with forage quality parameters than the original spectra or other spectra transformations (derivatives and continuum removal band depth). The normalisation process reduced the temporal noise of the data, which strengthened the relationship between spectral features and forage quality. Further, two correlograms from CP and ADF showed an inverse pattern, confirmed by the negative correlation between them. Nevertheless, both forage quality parameters had a stronger relationship in the red-edge and the near-infrared regions of the spectrum.
Linear regression models with spectral features (single waveband, NDSI, SR) obtained a maximum adj.R 2 of 0.42 for CP and 0.34 for ADF. The best correlated NDSIs comprised wavebands in the chlorophyll absorption features in the visible domain of the spectrum (470-518 nm and 550-750 nm), which relates to N and other biochemicals in green canopies [14]. These results paralleled with the results from previous studies that come from the same geographical region [15,18]. The study from Safari et al. [32] obtained linear regression models with NDSI with R 2 of 0.58 and 0.49 for CP and ADF, respectively. Moreover, the results from Biewer et al. [15] reported R 2 values of 0.33 and 0.13 R 2 for CP and ADF models with SR, respectively. Nonetheless, the results from the studies mentioned above were obtained with field spectroscopy data, which had a more comprehensive spectral range (305-1700 nm) than the camera utilised in the present study.
Several predictive modelling algorithms were tested to identify the best algorithms to estimate CP and ADF from the full spectral data, as no single consistent algorithm was shown to surpass all the given circumstances every time in a study by the authors of [61]. Moreover, model consistency was evaluated by training and testing 100 models using 100 different random train and test data sets, which allowed to disclose the model performance irrespective of the calibration data set. Except for the PLSR, other tested predictive algorithms (GPR, RFR, SVR, and CBR) proved promising for the underlying data. The SVR and CBR models for CP and ADF estimation showed the maximum precision (highest R 2 p ) and prediction accuracy (lowest nRMSE p ) followed by the RFR, GPR, and PLSR models, respectively. SVR and CBR predictive modelling algorithms were previously utilised to estimate water quality parameters based on spectral data [62,63], however, to our knowledge, such algorithms have not been employed so far to estimate forage quality parameters.
According to the literature, the PLSR and RFR were the most prominent predictive modelling algorithms in forage quality parameter estimation using spectral data. For example,  report on PLSR models, which obtained nRMSE values of 8.5% and 7.3% for CP and ADF, respectively. However, PLSR models resulted in the lowest accuracies both for CP and ADF in this study. Moreover, Pullanagari et al. [34] achieved an nRMSE of 11.2% for CP with the RFR model, and Singh et al. [33] obtained an nRMSE of 21.7% for ADF with the same modelling algorithm. It is noteworthy that the studies mentioned above only tested one predictive modelling algorithm; thus, no conclusions are possible considering the comparison with other algorithms.
Both CP and ADF estimation models resulted in less than 15% relative prediction error. However, CP estimation had a slightly lower relative error (median nRMSE p = 10.6%) than ADF estimation (median nRMSE p = 13.4%). A similar relative error pattern for CP and ADF estimation models was obtained in previous studies that utilised field spectroscopy data [16,18,64]. Moreover, the data points from the IMG were clustered out and acted as the driver of the CP model according to the plot of fit ( Figure 9a). Although, ADF data points in the observed against predicted plot did not highlight a similar pattern (Figure 9b). The high variation of CP between the cuts in IMG due to intensive management practice might be the reason for the mentioned pattern. However, the data points from other grasslands were mostly grouped in both plots because they almost experienced similar management practice.
Remarkably, the wavelength (718 nm) that obtained the highest linear correlation with CP turned out to be the most fundamental wavelength for the final CP and ADF models ( Figure A1). However, the identified first most important wavelengths were not related to any foliar chemical absorption feature, such as chlorophyll, protein, or cellulose [65]. According to Curran [65], the absorption features that relate to CP and ADF can be found at wavelengths in the shortwave infrared region of the spectrum (1400-3000 nm). Singh et al. [33] also confirmed that important wavelengths for ADF were in the shortwave infrared region of the spectrum. The camera used in the present study did not have the sensitivity to the shortwave infrared regions, even though the study showed a superior accuracy for grassland CP and ADF.
Imaging spectroscopy provides the advantage to generate maps by direct application of the model to recorded spectral images without the need for spatial interpolation. The prediction maps in our study displayed a substantial spatial variation in forage quality on the grasslands. While there is an inverse relationship between CP and ADF at the early cutting date of the MHML grassland, it vanishes with the increasing maturity of plants ( Figure 12). Thus, biomass from such grasslands (MHML, NSGL) may display a much more considerable variation in CP than in ADF content, which probably depends very much on the degree of Lupinus polyphyllus invasion. Such information may help farmers to identify areas with contrasting forage quality, which is a prerequisite for targeted site-specific management of larger fields, for example, to determine a sound scheduling of mowing or grazing activities. Further, biomass from grassland areas with contrasting quality may be collected and ensiled separately to match the requirements of animals, which usually differ sharply among different animal species and performance categories.
To summarise, predictive modelling algorithms allow an adequate forage quality estimation, regardless of the grassland type and cutting regimes applied. Ultimately, fine-tuning of the calibrated models with data from further diverse grassland types could increase the robustness of the models to generate forage quality maps for grasslands with different vegetation types and management practices.

Conclusions
The present study aimed to estimate CP and ADF concentrations of a multitude of grasslands with different vegetation composition and cutting regimes applied based on UAV-borne imaging spectroscopy data. It was demonstrated that the resulting models could accurately estimate CP and ADF irrespective of the grassland type, and that model accuracies are in the same range as those obtained with the use of field spectroscopy. With such models, forage quality maps can be generated, which may provide a benefit for practical grassland farming. Eventually, continued tweaking of the models with data from additional grasslands will be a way forward to ultimately obtain an accurate generalised model with a validity range that covers most grasslands within a broader region.