Understanding Growth Dynamics and Yield Prediction of Sorghum Using High Temporal Resolution UAV Imagery Time Series and Machine Learning

Abstract

Unmanned aerial vehicles (UAV) carrying multispectral cameras are increasingly being used for high-throughput phenotyping (HTP) of above-ground traits of crops to study genetic diversity, resource use efficiency and responses to abiotic or biotic stresses. There is significant unexplored potential for repeated data collection through a field season to reveal information on the rates of growth and provide predictions of the final yield. Generating such information early in the season would create opportunities for more efficient in-depth phenotyping and germplasm selection. This study tested the use of high-resolution time-series imagery (5 or 10 sampling dates) to understand the relationships between growth dynamics, temporal resolution and end-of-season above-ground biomass (AGB) in 869 diverse accessions of highly productive (mean AGB = 23.4 Mg/Ha), photoperiod sensitive sorghum. Canopy surface height (CSM), ground cover (GC), and five common spectral indices were considered as features of the crop phenotype. Spline curve fitting was used to integrate data from single flights into continuous time courses. Random Forest was used to predict end-of-season AGB from aerial imagery, and to identify the most informative variables driving predictions. Improved prediction of end-of-season AGB (RMSE reduction of 0.24 Mg/Ha) was achieved earlier in the growing season (10 to 20 days) by leveraging early- and mid-season measurement of the rate of change of geometric and spectral features. Early in the season, dynamic traits describing the rates of change of CSM and GC predicted end-of-season AGB best. Late in the season, CSM on a given date was the most influential predictor of end-of-season AGB. The power to predict end-of-season AGB was greatest at 50 days after planting, accounting for 63% of variance across this very diverse germplasm collection with modest error (RMSE 1.8 Mg/ha). End-of-season AGB could be predicted equally well when spline fitting was performed on data collected from five flights versus 10 flights over the growing season. This demonstrates a more valuable and efficient approach to using UAVs for HTP, while also proposing strategies to add further value.

1. Introduction

In the last 20 years, nucleic acid sequencing techniques have driven major advances in crop genomics, genetics and molecular biology [1,2]. Discovery science and crop improvement are now often limited by the speed and ease of obtaining large amounts of phenotypic information from crop trials [2,3,4]. Therefore, the development of methods for high-throughput phenotyping (HTP) of crops to reduce manual labor in the field and to better characterize phenotypes is increasingly required [3].

Recent technological advances provide new opportunities for the use of unmanned aerial vehicles (UAV) as a low-cost platform for carrying sensors that will deliver high spatial, temporal, and spectral resolution imagery to generate precise information about the interaction of solar radiation and vegetation [5]. In recent years, UAV-based structure-from-motion (SfM) techniques have been rapidly adopted to estimate traits such as canopy height [6,7] and yield [8]. Remote sensing from multispectral and hyperspectral sensors and image analysis techniques have been utilized to monitor nutrient status [9,10], above-ground biomass (AGB) [11,12], leaf area index [13,14], canopy cover [15], and senescence rate [15,16]. However, most analyses have been limited to one, or a small number, of sampling dates that are often focused towards the end of the growing season.

The tendency to focus on late-season data collection is logical given strong interest in economically important traits such as final yield or AGB. In addition, late-season phenotypes are the integrative result of all the developmental progressions and environmental conditions that occur throughout a growing season. However, there are numerous challenges associated with late-season phenotyping. Foremost, information on how the traits of interest vary across populations or treatments is often not available until very late in the season, or even after the season is complete if labor intensive destructive harvests are performed. This slows the rate of learning or discovery because material cannot be selected for more detailed analysis or breeding until the following growing season [2]. Second, damage to plants late in the growing season can introduce unwanted variation into phenotypic data. For example, while the heritability of plant height in a population of biomass sorghum increased over the first three months of the growing period, it declined in the fourth month when lodging occurred [17]. Consequently, it has been proposed that phenotyping secondary traits (e.g., green coverage, senescence rate, or crop height) that are correlated with the main target trait could avoid these problems [18].

Phenotyping as early as possible in the growing season allows more time to make decisions on subsequent sampling or selections for breeding prior to harvest for the next generation of testing [18]. Explicitly considering variation in crop phenotypes through time by repetitive monitoring can provide even greater quality data and insights about genetic and environmental variation [16,19,20,21,22,23,24]. However, collecting information repetitively over time requires a significant increase in effort and cost. So, evaluating the optimal intensity of sampling is of particular interest.

The high biomass yield potential, resource use efficiency and energy content of Sorghum bicolor (L.) Moench, along with the ability to grow on marginal lands, have made it a model C4 plant for bioenergy research worldwide [25,26,27,28]. However, sorghum performance could be further enhanced if large populations of diverse germplasm can be screened more easily using non-destructive and precise HTP methods [29]. Accurately assessing AGB is especially important and challenging for biomass crops such as photoperiod-sensitive sorghum [29]. Many genotypes achieve a final crop height of > 4 m [30]. Consequently, manual phenotyping is difficult and laborious, while mechanical harvesting requires specialized, expensive machinery [31]. With respect to HTP, the height of the crop also challenges the use of remote sensing approaches that use some ground vehicles [32,33]. In contrast, using UAVs to estimate AGB of sorghum is not hampered by the physical challenges of moving over or through the crop canopy [29]. Nevertheless, UAV-based structure from motion techniques and multispectral imagery for yield prediction have been applied to bioenergy crops much less than the major cereal crops [11,34]. Despite the rapid rate of growth achieved by these crops, most strategies have relied on single-date [29,34] or multi-date stacked imagery [5,35], meaning that information about the growth dynamics and temporal trajectory of the crop have not been exploited as a proxy of productivity. This brings opportunities to unanswered questions, e.g., what is the predictive value of repetitive screening of capturing plant growth dynamics, and how and when can we maximize information gain using UAV systems in the season?

In summary, the current work aims to investigate how SfM techniques and multispectral remote sensing data collected throughout crop development can be used to predict end-of-season AGB and to improve the efficiency and strategic use of UAVs as an HTP tool. Advances in this area are very broadly applicable because crop biomass is a core trait in many research and crop improvement programs. This includes the investigation of productivity, resource use efficiency, microbiota interactions and resilience to abiotic and biotic stresses. In the present study, the application is a characterization of genetic diversity in biomass production of sorghum as a biofuel feedstock. The nature of that variation and the high productivity of sorghum makes it a compelling test subject. Approaches to advancing the accuracy and efficiency of UAV-based phenotyping involved testing (1) the relative value of trait data from single dates versus traits describing dynamic growth as predictors of final AGB, (2) the influence of the temporal resolution of data collection (i.e., number of flights) on the prediction of final AGB, (3) how early in the season the prediction of final biomass could be achieved without a loss of accuracy, and (4) the relative value of geometric versus spectral trait data at different stages of crop development.

2. Materials and Methods

The overall workflow (Figure 1) involved: a) data collection on multiple dates from a large field trial of diverse biomass sorghum using a UAV carrying a multispectral (MSI) sensor, b) image preprocessing to produce orthophotos, (c) estimation of geometric and spectral features for each subplot on each sampling date, (d) spline curve fitting for each subplot, (e) extraction of geometric and spectral features from spline fits to produce time-point and dynamic growth variables for each subplot, (f) model training and cross-validation, (g) variable importance determination, and (h) model validation.

2.1. Experimental Design and Biomass Harvesting

The experiment used an augmented incomplete block design to grow 864 diverse accessions of biomass sorghum in single, four-row plots, plus six additional accessions in 16 four-row plots distributed across each of the 16 blocks, as has been previously described in detail [36,37]. The experiment was planted on May 31, 2019 at a site located at the University of Illinois Energy Farm research facility, Urbana-Champaign (40.065789°N, −88.208477°W; Figure 2a,b). Plots were harvested for AGB between September 14 and 17 using a four-row Kemper head attached to a John Deere 5830 tractor (Figure 2d). The total AGB and moisture content was determined using a plot sampler that had a near-infrared sensor (model 130S, RCI engineering). Plants in each plot were cut with the Kemper head and the harvested material entered the upper hopper of the sampler to be weighed. The material was then fed into a blower system and delivered through a spout where the near-infrared sensor is located to determine the moisture content based on established calibrations (Figure 2d). AGB in dry tons per hectare was calculated as follows: AGB dry tons/ha = total plot wet weight (kg) × (1-plot moisture) / (plot area in square meters/10,000). Strong late-season winds produced lodging damage at some specific locations in the field, making mechanical harvesting difficult. For this reason, data from the 782 undamaged plots were used for the analysis.

2.2. UAV System, Data Collection, and Pre-Processing

A multispectral camera (Rededge-M, Micasense, Seattle, WA, USA) was used for imaging to acquire both photogrammetric and multispectral traits. The camera is a global shutter sensor with five spectral bands in the blue (465 to 485 nm), green (550 to 570 nm), red (663 to 673 nm), rededge (712 to 722 nm), and near-infrared (820 to 860 nm) regions of the electromagnetic spectrum. The fast shutter speed of the sensor reduced blur and artifact effects in the images, which resulted in high-quality 3D reconstruction of the canopy through the season. Intrinsic optical parameters and spectral calibration correction was implemented in Metashape software (Agisoft, St. Petersburg, Russia, 2016). A standard Micasense calibration panel was imaged on the ground before and after each flight for spectral calibration. Known albedo values of the panel for each spectral band were utilized in an empirical calibration procedure as follows:

R_sunlit(i) = (S_sunlit (i) / S_ref (i)) ∗ R_ref(i)

(1)

In Equation (1), S_sunlit and S_ref are radiance values from sunlit leaves and the reference Micasense panel, respectively, R_ref refers to the reflectance of the calibration panel (calibrated, known value) provided by Micasense, and R_sunlit is the absolute reflectance of sunlit leaves in band i.

The aerial platform utilized was a Matrice 600 Pro hexacopter (DJI, Shenzhen, China), equipped with a Gremsy T1 gimbal (Gremsy, Ho Chi Minh, Vietnam) used to mount the multispectral sensor (MSI), shown in Figure 2c. Flights were conducted 10 times in the season between 21 and 98 days after planting (DAP; Figure 3a) under clear sky conditions, and between 11:00 and 13:00 local time. Flight planning allowed for 90% forward and 80% side overlapping in image acquisition following the Metashape software (Agisoft, St. Petersburg, Russia, 2016) recommendations for successful 3D reconstruction, flying at an altitude of 40 m above ground level, resulting in a ground sampling distance (GSD) of 2.5 cm/pixel. A Trimble R8 global navigation satellite system (GNSS) integrated with CORS-ILUC station was utilized for RTK (real-time kinematic) surveying of nine ground control points (GCPs) placed along the edges and the center of the experiment, with a distance of 60 to 65 m between each other. All this information was utilized for precise geo-referencing and co-registration of orthophotos between dates of data collection.

The excessive greenness index (EXG) [38] was utilized to binarize vegetation and background pixels to prevent contamination and misrepresentation of vegetation in the images. This step is necessary to mask out pixels with abnormal spectral values and regions of the image including shaded leaves and background. Pixels classified as background were masked out and not considered for the extraction of summary statistics from each plot. Green coverage (GC) is defined as the ratio of green versus total pixels in each plot. The earliest flight of the season on the day after planting (DAP) 21 was utilized as the ground level reference or digital terrain model (DTM) for the calculation of the absolute canopy height or crop surface model (CSM) on subsequent sampling dates. The results were consistent between early-season DTM and CSMs through the season. Plot polygon boundaries were delineated in QGIS v3.14.1-Pi (QGIS Geographic Information System; Open-Source Geospatial Foundation Project). Zonal statistics for each feature (

Table 1. Description of UAV-based geometric and spectral features utilized, their corresponding equations, and case studies reporting them.

Features Variables	Description	Formula	Reported by
CSM	Geometric	CSM = DSM − DTM	[39]
GC	Geometric	(n pixels green/n pixels total) ∗ 100	[29]
NDVI	Spectral	(NIR − Red/NIR + Red)	[40]
NDRE	Spectral	NIR − Rededge/NIR + Rededge)	[10]
WDRVI	Spectral	(0.1 ∗ NIR − Red/0.1 ∗ NIR + Red)	[41]
NGBDI	Spectral	(Green − Red/Green + Red)	[42]
EXG	Spectral	(2 ∗ Green − Red − Blue)	[38]

) were calculated at each plot, with the 90% trimmed mean value of each feature being utilized in further analytical steps. The resulting estimates of canopy height from UAV-based images were validated against manual measurements at DAP 72 and showed a high degree of correspondence (R² = 0.84) (Figure S1).

2.3. Temporal Integration via Spline Fitting

Spline fitting enabled the integration across time of information from multiple single flights, while also correcting and harmonizing for the uneven timing between flights that was caused by inclement weather conditions (e.g., precipitation or cloud cover). Consequently, an evaluation of data from individual time-point versus dynamic features of crop growth throughout the entire growing season was made possible.

Due to its particular photoperiod sensitivity [43], biomass sorghum continues to produce new vegetative biomass throughout the growing season at this location. Therefore, the growth trajectory of the crop does not clearly plateau as annual crops with a strong mid-season vegetative-reproductive transition do. This continued vegetative growth makes the traditional approach of fitting Gompertz or logistic curves to data challenging for this crop. Cubic smoothing splines have been successfully utilized [44] to characterize the growth of plants over time. It is a flexible curve-fitting technique that integrates cubic splines and curvature minimization to create an effective data modeling tool that preserves the relevant signal but removes uninformative noise in the temporal profile [45]. Specifically, a penalized criterion known as the penalized sum of squares (PSS) is defined for cubic smoothing splines (Equation (2)). The combination of both functional constraints makes smoothing splines a suitable and flexible alternative to characterize the growth dynamics of biomass sorghum.

$$\mathrm{PSS}={\displaystyle \sum }_{i=1}^n{\left(yi-\mathrm{S}\left(xi\right)\right)}^2+\lambda {\displaystyle \int }{\mathrm{S}}^{″}\left(x\right){)}^{2}dx$$

(2)

In Equation (2), the penalized sum squares criterion is defined, the integration occurs over the range of x and λ is a tuning parameter. The first portion of the equation defines a first functional constraint that minimizes the squared error between the data and spline; the second portion of the equation considers a second constraint that directly penalizes the curve flexibility.

The smooth.spline function of the stats library in R software (R Core Team, Vienna, Austria, 2019) was utilized to implement spline fitting for each plot. The function used the input data as knots and used generalized cross-validation (GCV) to prevent under- or overfitting of the spline lambda parameter [46,47] and find an optimal level of smoothing. The minimum number of knots used was four, which is coincident with the minimum number of flights considered for the spectral features in flight pattern b (Figure S3b). The function was set to return an approximation of a continuous solution on a daily time-step over the interval of the flights (e.g., Figure 4a,b, Figures S2a,b and S3a,b) as well as the first derivative of that solution (e.g., Figure 4c,d, Figures S2c,d and S3c,d) (Table 1). Finally, the solution and its first derivative were returned at eight specific predefined dates, as described in Figure 5. In order to test the effect of the sampling frequency, the spline fitting procedure described above was applied to data from five dates (Figure 3b), as well as the full set of 10 dates (Figure 3a).

After spline fitting, the value of each of the seven features (Table 1) and the first derivative for each trait was extracted at eight predefined DAPs (30, 40, 50, 60, 70, 80, 90, 98) for each subplot (e.g., Figure 5). With eight dates and seven features (Table 1) assessed as both single time-point and dynamic traits, 112 traits were generated in total. Traits were named using the convention of “type of feature_DAP_type of variable”. For example, “NDRE_40_tp” describes NDRE on DAP 40 as a single time-point trait, while “DSM_30_slp” describes the first derivative of canopy height at DAP 30 as a dynamic growth trait.

2.4. Machine Learning Algorithm and Variable Importance

Random Forest (RF) was selected over other machine learning methods because, as a non-parametric tool, it (1) does not assume the data is normally distributed and (2) does not make assumptions about the form of associations between predictor variables and the response variable, while (3) simultaneously using multiple predictors [48]. Furthermore, as an ensemble of trees, RF is more capable in terms of prediction power and reduces overfitting compared to simple tree learners [49]. RF was used for two purposes: (i) to assess the relevance of features on prediction, and (ii) as a tool to predict AGB.

Conditional inference trees were used to build the ensembles (forests) using the cforest function [50] via caret library [51] in R software. The use of permutation importance via cforest function ensures feature selection was not biased towards predictors with many possible splits [48]. This variable importance procedure [49] has been demonstrated to reduce bias compared with other alternatives [52]. Furthermore, since the dataset includes correlated features (Figures S4 and S5), variable importance was evaluated using a conditional permutation test to minimize the overestimation on importance scores of correlated features [48].

2.5. Prediction Models

Three model configurations were tested for the prediction of end-of-season AGB {y} from a series of input variables {x₁, x₂, …x_n} at each predefined DAP, with the aim of quantifying the relative contribution and AGB predictability of time-point and dynamic features throughout the season:

Model 1:

(time-point feature x₁+…+ x_k) = {AGB}

Model 2:

(dynamic feature x₁ +…+ x_n) = {AGB}

Model 3:

(time-point feature x₁ +…+ x_k + dynamic feature x₁ +…+ x_n) = {AGB}

Data from the 782 subplots were randomly assigned for either training (70%) or testing (30%) the model. A tenfold repeated cross-validation resampling procedure was utilized for model training. During the training process, parameters mtry (number of variables at each tree node) and ntree (number of trees) were optimized using a grid search approach, which found the optimal parameter combination by minimizing the root mean square error (RMSE). The out-of-bag sample was not used to generate the actual tree, but was used to determine the permutation importance score of each variable by accounting for the increase in mean square error of the prediction once a variable was removed from the model [49]. Each model fitting was iterated 10 times using a random training and testing partition as a measure of uncertainty of the model at each predefined DAP. At each iteration, the best-tuned model was exposed to the unknown test dataset to evaluate how the trained model generalized the prediction of AGB.

The root mean square error (RMSE), mean absolute error (MAE), absolute error (AE), and coefficient of determination (R²) were utilized as evaluation metrics to determine the performance of the models. Each metric is described in Equations (3)–(6):

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{\left(yi-\widehat{y}i\right)}^2}$$

(3)

$$\mathrm{MAE}=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n\left|yi-\widehat{y}i\right|$$

(4)

$$\mathrm{AE}=100*\frac{{{\displaystyle \sum }}_{i=1}^n\left|\left(yi-\widehat{y}i\right)/yi\right|}{n}$$

(5)

$${\mathrm{R}}^2=1-\left[{{\displaystyle \sum }}_{i=1}^n\left(yi-\widehat{y}i\right){}^2/{{\displaystyle \sum }}_{i=1}^n(yi-{\overline{yi}}^2)\right]$$

(6)

In Equations (3)–(6), $y$, $\widehat{y}$, and $\overline{y}$ are the observed, predicted, and observed mean AGB values of the ith plot, and n is the total number of samples within the study area.

3. Results

3.1. Variable Importance

Analysis of variable importance focused on model 3 where all variables were available to inform predictions of end-of-season AGB. The relative importance of different variables to predictions of end-of-season AGB by model 3 were very similar for datasets produced from 10 flights or five flights distributed across the growing season (Figure 6). Geometric features describing the height (CSM) and ground cover (GC) of the crop canopy were generally more influential in predictions of end-of-season AGB by model 3 than spectral features such as EXG, NDRE or NDVI. This was the case regardless of which timepoint of data predictions were based on, and for both single time-point (*_tp) traits or the rate of change of the trait value on a given date (*_slp). However, the importance of individual features varied over the growing season (Figure 6). In general, GC was more important prior to canopy closure, which occurred around 60 DAP (Figure S2); dynamic variables describing the rate of change of traits were more important early in the season rather than late in the season (Figure 6).

Specifically, on the earliest date (30 DAP) the rate of height growth (CSM_slp), the rate of gain of ground cover (GC_slp), and the absolute GC on that date (GC_tp) were the most important variables for prediction of AGB by model 3. Over the subsequent month, the canopy height on given date (CSM_tp) became progressively more important. Meanwhile, CSM_slp and, to a lesser degree, GC_tp and GC_slp became progressively less important. Late in the season, from 70 to 98 DAP, CSM_tp was the only very important variable. The spectral features of NDRE_slp, NDRE_tp, EXG_slp and EXG_tp were of only moderate importance throughout the growing season, but did explain residual variance after geometric features were accounted for (Figure 6).

3.2. AGB Prediction

Using data from five rather than 10 flights did not significantly reduce the performance of any of the three categories of models or alter their relative performance in terms of RMSE, MAE, and R² (Figure 7). Moreover, predictions based on data from both flight patterns evidenced a similar prediction stability, as indicated by their standard deviations (Figure 7d–f) and for predictions based on data from 10 flights (Figure 7a–c).

Combining single-time-point and dynamic variables in model 3 generally resulted in the best and most consistent performance in terms of RMSE, MAE and R² throughout the growing season (Figure 7). Error was lowest (RMSE = 1.7–1.8 Mg/ha; MAE = 1.3–1.5 Mg/ha) and the proportion of variance explained was greatest (R² = 0.59–0.63) for model 3 between DAPs 40 and 50, i.e., shortly prior to canopy closure. Model 1 performed essentially equivalently to model 3 when based on single-time-point data collected after canopy closure (>70 DAP). However, predictions of end-of-season AGB by model 1 were less accurate and explained less variance than model 3 when based on early-season data. In contrast, predictions by model 2 were better when based on dynamic variables measured early versus late in the growing season. However, even model 2 did not perform as well as model 3 prior to canopy closure.

3.3. AGB Prediction for Six Accessions Tested in Highly Replicated (n = 16) Plots

Mean end-of-season AGB varied from 21 to 27.6 Mg/ha across six accessions grown in highly replicated (n = 16) plots (Figure S9). However, equivalent or greater variation in AGB was observed among the 16 replicate plots of each accession (Figure S9). Consistent with results for the population as a whole (Figure 7), combining single time-point variables and dynamic variables in model 3 produced accurate and consistent results at both the early season stage (DAP 40) and after canopy closure (DAP 80) (Figure 8). By contrast, and more clearly marked in the reduced number of flights pattern, single time-point variables model 1 produced higher mean AE early in the season (Figure 8a,c) in three and four of the six genotypes of each flight patterns, respectively. Aligned with the whole population prediction, using only dynamic variables in model 2 produced higher AE late in the season (Figure 8), where the predicted AGB does not align accordingly with the observed values (Figure S9b,d).

Even with model 3, there is a tendency towards underestimating the end-of-season AGB of the most productive genotypes (e.g., PI148089; Figure S9b,d). However, that was less the case early in the season (DAP 40) than late in the season (DAP 80). A general pattern of bias leading to underestimations of AGB in high-yielding genotypes was evident in the population as a whole (Figure S8). The variance among replicate plots of a given genotype was often substantially less for model predictions of end-of-season AGB than it was for data from the mechanical harvest (Figure S9). Again, a reduction in the number of flights from 10 to five did not substantially affect how models 1, 2, or 3 performed in predicting AGB of these highly replicated accessions, as seen in Figure S9.

3.4. Trait Relationships

Patterns of pairwise correlations among traits were very similar for datasets based on spline fitting of data from 10 flights per season (Figure S4) compared to five flights per season (Figure S5). So, interpretation of the results will focus on the use of the full dataset from 10 flights. Trait relationships varied significantly over the course of the growing season (Figure S4). Early in the season (DAP 30, 40, 50, 60), many measures of growth rate and absolute size were associated with each other; i.e., there were positive correlations among the single time-point spectral traits (NDVI_tp, NDRE_tp, WDRVI_tp, NGRDI_tp), the single time-point geometric traits (CSM_tp, GC_tp), the rate of canopy growth both vertically (CSM_slp) and horizontally (GC_slp), as well as end-of-season AGB (Figure S4). In contrast, on DAP 30 and DAP 40, the initial rates of gain of spectral indices were correlated with each other, but not with the single-time-point spectral or geometric traits describing absolute plant sizes. Additionally, during the period when canopy closure occurred for different genotypes (DAP 50, 60, 70), the rates of gain of ground cover (GC_slp) and spectral traits (NDVI_slp, NDRE_ slp, WDRVI_ slp, NGRDI_ slp, EXG_slp) became negatively correlated with single time-point measures of geometric and spectral variables (CSM_tp, GC_tp, NDVI_tp, NDRE_tp, WDRVI_tp, NGRDI_tp; Figure S4). Late in the growing season, there were no significant correlations among the single time-point traits and dynamic traits describing the rate of change of either geometric or spectral features. However, there was still a significant correlation between canopy height and end-of-season AGB.

4. Discussion

This study successfully addressed its goals by demonstrating that early- and mid-season measurement of dynamic growth traits by UAV can be combined with individual time-point trait data to facilitate improved prediction of end-of-season AGB (Figure 6 and Figure 7). Specifically, data on the rate of change of geometric and spectral features extracted from multispectral imagery of biomass sorghum canopies resulted in more accurate predictions of end-of-season AGB at much earlier dates in the growing season. In addition, end-of-season AGB could be predicted equally well when spline fitting was performed on data collected from five flights versus 10 flights over the growing season. The resulting analysis revealed limitations and opportunities for HTP of genetic variation in the productivity of photoperiod sensitive biomass sorghum accessions, which presents special challenges as a result of its delayed transition to reproductive development and high AGB.

Geometric and spectral information have proved to be relevant descriptors of biomass in annual grain crops such as rice [53], barley [54,55], and wheat [56]. Spectral information has proved to be an important proxy for AGB prediction in annual grain crops in these studies, but prediction accuracy significantly increased when integrating both types of traits in the same model [54,56]. It was also pointed out that geometric features provided more robust AGB prediction results in barley [55] and wheat [56] or provided both similar prediction power in rice [53]. In our study, geometric information is consistently the more influential descriptor of final biomass than spectral information, regardless of whether it is being considered as a trait reported at a single time-point or as a dynamic trait describing a rate of change over time (Figure 6). This is consistent with canopy height previously being identified as a trait that displays significant intraspecific variation, with high heritability, in other species [53,54,55,56,57,58]. This is probably, in part, because it does not suffer from signal saturation as spectral information does late in the growing season. Biomass sorghum shows a rapid increase in ground cover and vertical growth, such that canopy closure and heights of ~1.5 to 2.5 m are achieved at around 50 to 70 DAP (Figure 4 and Figure S2). The power to predict end-of-season AGB was greatest around DAP 50, presumably because this gave the genetic variation present in the large diverse population of accessions the greatest period of time to express itself in AGB before occlusion associated with canopy closure made it harder to image the size and structure of vegetative structures. Demonstrating the value of measuring early-season dynamic variables is a departure from the conclusions of previous studies where only single time-point data were analyzed and the best predictions of biomass were found using mid- and late-season data [18,32]. However, there is still the common finding that single time-point height data is a key predictor of end-of-season AGB when measured late in the growing season. Moreover, it is important to recognize that stressful or extreme weather conditions late in the season could cause end-of-season AGB to deviate from what is predicted from early-season indicators of productivity. So, additional analyses across multiple years or locations will be needed to establish how robust the findings of this study are across environmental variation.

Highly replicated plots (n = 16) of six genotypes provided a valuable test of HTP methods that complemented the broader survey of the full population of biomass sorghum. The results were very similar to what is described above in terms of the value of geometric versus spectral traits early versus late in the growing season and with data from 10 flights versus five flights. It was notable that the underestimation of end-of-season AGB in the highest productivity line by model 3 (time-point + dynamic variables) was less in early-season predictions based on dynamic traits than in late-season predictions driven mainly by single time-point height data (Figure 8). It could be argued that the general underestimation of AGB in high-productivity accessions is directly related to the distribution of the AGB data with regard to extreme values [5]. There was substantial variability in end-of-season AGB across the 16 replicate plots of these six accessions. Some of this will have resulted from spatial variability across the field, but there is also likely significant variation associated with inconsistency in biomass estimates from the forage harvester. Variability in predictions of AGB from the UAV data were notably less variable across replicate plots. Additional experiments will be needed to demonstrate whether this is a consequence of the RF modeling approach or results from UAV measurements leading to more reproducible AGB estimates than destructive harvest.

The potential to roughly predict variation in end-of-season AGB from dynamic traits at mid-season or earlier could be very valuable for a number of reasons. A subset of genotypes with desirable characteristics could be selected for: (1) seed production by crossing or selfing well in advance of key reproductive events; or (2) additional phenotyping by manual methods or other HTP techniques. In the case of short-day photoperiod sensitive sorghum, plants might even be transplanted out of the field into controlled environment conditions, where flowering could be induced earlier than would occur naturally in a mid-latitude location such as the U.S. Midwest [59]. In both cases, labor and expense could be saved within the field season and some activities could be performed without the usual need to wait for an additional generation of plants to be grown. The value of accelerating breeding cycles is well recognized in many major food crops [60], and the improvement of bioenergy crops would also benefit from this strategy.

Many previous studies indicate that spectral vegetation indices lose sensitivity when biomass accumulation is high and canopies are dense [41,42]. This is consistent with the negative correlations that were observed between the first derivative of spectral features and measures of absolute canopy size at mid-season (Figure S4). The limited value of spectral features to predictions of end-of-season AGB in the current study probably also resulted from limiting the trial to the study of genetic variation in a single environment. This is a key context in which HTP needs to be tested, but where there is not as much phenotypic variation to exploit as in experiments where environmental or management treatments such as variable nitrogen fertilization [29] or drought stress [61] are also imposed. An important practical consequence of being able to predict end-of-season AGB from dynamic geometric features early in the season is that it indicates that low-cost RGB cameras can be utilized as lower-cost alternatives to multispectral sensors without significant compromises on the value of the data collected [62].

Reducing the number of data collection flights, from 10 to five dates spaced across the growing season, did not penalize the sensitivity and ranking of variables’ importance through the season. Since the crop reaches full green coverage in a very short period (Figure S2) flights earlier than DAP 50 become critical to differentiate final productivity between genotypes. This highlights that a reduction of the number of flights can be implemented, but timing of the flights is key for specific variables. To overcome this challenge, prior understanding of the seasonal predictability of the trait of interest, as reported in this work, becomes critical for efficient use of the UAV in the season. By accounting for that, researchers will be able to prioritize collection of aerial data at stages of the crop with high information value, rather than collecting aerial data with marginal informative value, or by reconstructing the entire temporal season profile of the crop, as presented in this work.

While the method presented here provides valuable predictions of end-of-season AGB from early-season phenotype data, there is still significant unexplained variation in AGB. Sorghum accessions vary significantly in the number and diameter of stems [63,64]. Variation in those traits will be hard to detect in geometric or spectral features of aerial imagery after canopy closure. Therefore, it seems likely that aerial and ground-based imaging using rovers or push-cart mounted sensors [65,66,67] will need to be paired together to capture more of the traits needed to explain the full diversity of AGB in biomass sorghum.

5. Conclusions

This study demonstrated the use of high temporal resolution UAV imagery to understand the relative importance of dynamic and single-date static information throughout the season to predictions of final harvestable AGB. Geometric features were consistently the most informative variables of end-of-season AGB in the season. In particular, rapid rates of growth early in the season were positively associated with final productivity. In addition, by being able to identify stages of the crop with higher information gain, a more precise and efficient strategic use of the UAV as HTP tool in the season can be implemented. Spectral information was not a main descriptor of end-of-season AGB, which shows RGB sensors are expected to reach similar performance as multispectral sensors but at a lower budget. The number of visits to the field, which are an important component of the cost of operations of these platforms, can be significantly reduced with no information lost if data collection is timed appropriately. Future directions should explore the impact of the spatial resolution of imagery, and the complementary integration of the capabilities of UAV to cover large areas and high sensitivity of ground vehicles to characterize canopy architectural traits to overcome limitations associated with the occlusion of stems and the saturation of spectral features after canopy closure.