The Second Derivative of the NDVI Time Series as an Estimator of Fresh Biomass: A Case Study of Eight Forage Associations Monitored via UAS

In the current climate change scenario, the prediction and estimation of crop production in terms of biomass or yield are paramount tasks for the agricultural scientific community; these have become urgent and highly important objectives in the response on the global scale to the food demands of a growing world population [1].

Crops whose productivity depends on the amount of plant biomass usually belong among forage cultivars, which are considered the main feed source for ruminants. These animals are capable of transforming an initially poor-quality foodstuff, generally cereal-based forages with a low protein content, into high-quality products for human consumption, such as milk and meat [2]. Thus, the improvement of forage quality can be translated into the optimization of animal productivity; thereby, more inefficient and unsustainable practices, such as purchasing protein supplements, can be avoided [3]. In this context, intercropping systems, i.e., the simultaneous growth of two or more species [4] in which cereal and legume species are combined, have emerged as a sustainable alternative that may be used to increase the quality and quantity of forages [5]. In contrast to grain crops, these forages are usually harvested with a certain moisture content and have the best nutritional properties; thus, their yields are quantified primarily in terms of fresh plant biomass.

Among the many methods for producing crop statistics described by the Food and Agriculture Organization (FAO) [6], the yield and biomass statistics are usually gathered using crop cuts and/or farmer declarations. During highly expensive surveys, experts may sample subplots within the plot and measure the production by area. This may include visual estimations and the giving of questionnaires to farmers. Alternatively, farmers may be requested to provide a post-harvest estimation of the production in a given area, although visits to parcels and granaries are common. The complexity of estimating biomass over mixed crops is even higher, since intercrops involve two or more crops growing in the same field. In this context, remote sensing and modeling approaches are emerging as alternative procedures to crop production or yield [6].

The FAO has also promoted global initiatives to facilitate the study and applications of the relationship between crop yield and water use by publishing comprehensive guides about the generic topic of “irrigation and drainage”, namely, paper 24: “Crop water requirements” [7], paper 33: “Yield response to water” [8], and paper 66: “Crop yield response to water” [9]. All of these publications attempt to assist farmers and agricultural managers in implementing effective agricultural practices to enhance yield while preserving water consumption. In particular, the FAO paper 66 relies on the AquaCrop model to simulate biomass and yield [10]. AquaCrop belongs to a model group that addresses crop biomass productivity in relation to water availability [11]. Therefore, it is a so-called water-driven model based on the concept of water productivity (WP). The conceptual equation at the core of AquaCrop states that biomass production is proportional to the cumulative amount of water transpired [12]. Since its launch, AquaCrop has been at the core of many studies [13,14,15], including a specific version supported by geographical information systems and remote sensing [16,17]. In fact, the new developments are mainly oriented toward remote sensing data assimilation, which provides the missing spatial information required [18].

However, other larger groups of models exist, either carbon-driven or solar-driven, in which the growth (and, therefore, biomass) estimations are based on carbon or light assimilation, respectively. A ten-year review of AquaCrop [18] performed a thorough and up-to-date revision of different families of models, with many examples. Far from physical modeling, other statistical approaches to predict yield and water production are also nurtured by satellite imagery. Simple or multiple regression models relating remote sensing data features with productivity are still in use [19,20] since they are the simplest and easiest methods to compute, although the results are often inconsistent and not easy to generalize [21]. More sophisticated artificial intelligence techniques from the fields of machine and deep learning, which apply algorithms based on convolutional neural networks, support vector machines, random forest, etc., have been progressively implemented owing to the increasing availability of large and high-quality datasets [21]. The work of van Klompenburg et al. (2020) [22] presents a systematic literature review of these crop yield prediction alternatives.

As ground-based phenological observations are limited, phenology derived from remote sensing can be used as an alternative to parameterize phenological models [23]. Remote sensing data offer many advantages in crop prediction. First and foremost, the images provide a wide spatial range and scalability; they are spatially seamless and may fill in situ data gaps [24]. Second, satellite data can provide a synoptic overview of actual growing conditions and can be used to diagnose discrepancies from normal conditions [25]. Finally, remote data assimilation allows different alternatives: direct substitution of the model parameters, calibration/initialization, or sequential assimilation of algorithms and models [26,27]. In particular, in AquaCrop, the remote sensing inputs are usually assimilated as either indicators or surrogates of the model parameters and include the data (or their derivative products) of remotely sensed temperature, vegetation indices, leaf area indices, soil moisture, the fraction of photosynthetically active radiation, and many others [27,28,29,30,31]. Other examples use satellite data to calibrate AquaCrop inputs [28,32,33], and in other cases, the remotely sensed time series provides the temporal metrics of the growing cycle, such as the start, end, or length of the growing season [16,30].

Among the many remote sensing products, the normalized difference vegetation index (NDVI) [34] seems to be the most popular input in crop production simulation models. The NDVI is simple and easy to interpret and is readily available from most satellite providers [32]. In addition, all multispectral commercial cameras on board unmanned aerial systems (UAS) include red and infrared bands, which are the basis of its calculation [35]. Many examples of yield estimation through UAS observations have also been proposed [35,36,37]. These platforms provide a superhigh spatial resolution but present the disadvantage of their lack of an automatic revisit, as in the case of satellites. Occasionally, this might hamper the monitoring of the complete growing cycle, including the key moments of crop development or senescence.

All applications of the NDVI for yield, biomass, or water productivity are related to the well-known fit between the NDVI time evolution and the growing cycle and phenology of many crops [35,38,39]. Therefore, the NDVI has been used as an effective indicator of crop yield or plant biomass [20,40,41] from different perspectives, such as the direct correlation or the aforementioned assimilation modeling. The shape of the NDVI curve and the particular features of it, such as the integrated, maximum, and moving average or the relative range, have been used as synoptic indicators of biomass or production [19,20]. As another indicator, the second derivative of the NDVI curve has been used for estimating phenological information such as the start of the growing season [42], while its maximum has been related to the beginning of the green-up phase [23,43]. The inflection points resulting from the second derivative of the NDVI curve have inspired our hypothesis on the determination of the period in which the biomass is produced, together with the NDVI value itself.

The integration of NDVI values to estimate plant biomass, net primary production, or grain yield is one of the most popular approaches for the remote sensing community [54,55]. The main difficulty of the retrieval based on the NDVI integration lies in (1) choosing the precise time-critical crop growth stages that lead to the final biomass accumulation and (2) searching for these key phenological stages and identifying them over the NDVI time series within which the integration is to be performed. Regarding the first issue, the proposal was made to begin biomass accumulation during the late flowering and maturity stages, as was also suggested by Hassan et al. (2019) [35], rather than at the green-up of the initial development phases, as proposed by many other authors [30,40,54] who sought a sudden increase that might signal the onset of significant photosynthetic activity [55]. This selection is in line with Calera et al. (2004) [51], who selected the NDVI plateau stage coinciding with the linear growth phase as an estimator of the potential rate of matter accumulation. In addition, this choice fits the FAO66 condition of using Kc,Tr when the canopy cover is full, since the NDVI second derivative limits ensure a CC close to 100%.

Regarding the second question, i.e., how to identify the NDVI critical growth stages, thresholds based on the NDVI maxima, averages, or moving averages have frequently been proposed [56,57]. Although these thresholds possess the advantage of being easy to recognize, they are not systematic and can be affected by local conditions, such as vegetation, soil, and illumination [55], as well as by sensor biases. The use of NDVI derivatives, although less frequent in the literature, may overcome these problems. The first, second, and third derivatives were applied to determine the start of the season as the date of the maximum increase in the respective NDVI derivative curve [23,42]. In particular, the maximum of the second derivative has been related to the beginning of the green-up phase [43] or to the time when the majority of pixels are turning green [42]. Our data confirmed that the maximum of the second derivative indicates the downward concavity of the NDVI curve (end of March in Figure 3a,b) and thus also the onset of the green-up. However, our perspective is rather different. The biomass accumulation was produced during a later stage, coinciding with the local minimum of the second derivative and the maximum convexity of the NDVI curve (Figure 3). In particular, in our study, this point corresponded to 14 April. Our proposal agrees with Labus et al. (2002) [58] and Doraiswamy and Cook (1995) [59] in that the early-season NDVI parameters were not consistent indicators of the wheat yields, and the NDVI growth profiles showed a stronger relationship with the yield later in the season during the grain-filling stage. Conversely, at the end of the curve, the green biomass production decayed before the offset of the cycle, in the middle of May, coinciding with the latter minimum of the second derivative. This point would be readily identifiable at the end of the cycle and would allow farmers to harvest during an optimal time of production and forage conditions. The method may also provide an early warning about a potential harvest decline due to adverse weather or crop conditions. In addition, this early-season estimation could not only reduce resource input and environmental pollution, but also increase crop yield and the subsequent profits [60,61], as well as determine inputs such as nutrients, pesticides, and water in order to optimize the yield potential.

While the NDVI has repeatedly been used in AquaCrop as a surrogate of some parameters, including CC or crop coefficients [16,30,31,32], our approach considered the NDVI as a single synoptic indicator of crop vigor, biomass production, and plant status. Supporting this idea, it was shown that the NDVI may include different stresses, such as the impacts of fire, frost, or drought, during sensitive crop stages [25]. Therefore, following the water productivity fundamentals, the integrated NDVI was multiplied by WP* as a constant. In this first attempt, a value of

$$

W


P


*


=
18
 
g
r
/


m


2

(180 kg/ha) was chosen based on the literature. However, this value could be adapted to the specific characteristics of the forage associations. As a given example, taking into account that the NDVI maximum is 1, the expected maximum daily productivity would be 180 kg/ha. Furthermore, considering that the number of days ranged between 34 and 50 (depending on the association and the scale), the possible maximum fresh biomass production ranged from 6300 to 9000 kg/ha, which are considered feasible values of forage yield in a semiarid rangeland ecosystem such as the dehesa [62,63]. The actual results for the forage associations ranged from 4200 to 7400 kg/ha (Table 5) and were in line with the observed values and the different behavior observed with respect to the vetch-based and pea-based associations. The accuracy of these estimations in terms of the RMSD (excluding the failed date at the end of May) ranged from 700 to 2000 kg/ha, and the MAB ranged from 600 to 1900 kg/ha (all the values were positive), which represented error percentages ranging from 10 to 30% (Table 6 and Table 7). Although rather high, these numbers are in accordance with other crop biomass/yield estimations that used remote sensing and more complex approaches or models. Ajith et al. (2017) [64] used an integrated NDVI for rice and found error percentages ranging from 11 to 21%. Benedetti and Rossini (1993) [65] found an error from 10 to 19% using cumulative NDVI profiles in a regression model for wheat, and in a similar approach, Doraiswamy and Cook (1995) [59] found much higher and variable errors at the county/regional levels in two regions of the USA. Using the simple algorithm for the yield estimates (SAFY) model for maize with high spatial and temporal resolution remote sensing data, Battude et al. (2016) [66] found a relative error for dry biomass of 14%. Using the same model, Claverie et al. (2012) [67] found a relative error of 25% for maize and of 39% for sunflower. Ji et al. (2022) [68] used three machine learning algorithms for faba beans, obtaining yield estimation errors from 18 to 31%. Many other examples may be cited [69], although most of them present the error in terms of absolute errors of biomass (kg or tons per ha) for a given crop; so, they are hardly comparable to the forage associations studied here, which, in addition, are not frequently studied in the literature. In any case, the estimation may supply a rough estimate of forage biomass for livestock feeding that could otherwise have been difficult to draw from traditional surveys, as also stated by Benedetti and Rossini (1993) [65].

In addition to its simplicity, another relevant advantage that should be emphasized with regard to this approach is its objectivity. Occasionally, the parameterization of a model is applied in a relative manner, through comparison with other users or the same user over time, rather than applied as an absolute value [69]. This objectivity is sustained by the NDVI’s capacity to determine the status and vigor of the plant, which ultimately implies its biomass production capacity. While a higher NDVI value is associated with a faster growth rate and higher biomass accumulation during the vegetative stage [35], this higher NDVI is not always related to a higher grain yield. Nevertheless, there is a deep discussion in the scientific literature about the predictive capacity of the NDVI for biomass or grain yield depending on the developmental stage [40,65,70]. The approach performed reasonably well for fresh or green biomass, but this was not the case for the final biomass estimated at the end of May, when senescence had already initiated. The relationship between the biomass measured on 29 May and the NDVI on the same date for the eight associations showed good agreement in terms of R (Table 4), but substantial differences between the observed and the estimated weights were found on this date (Table 2 and Table 5, respectively). Hence, the MAB and RMSD were unacceptable. It should be highlighted that the focus of the work was the estimation of green forage for livestock feeding at the time of the highest nutritional quality of the forage crops. Hence, the last date was beyond the focus of the research.

In addition, the objectivity of the approach was built upon a new, systematic way to determine the key stages through the second derivative, although much more research should be conducted with other crops, different edaphoclimatic conditions, and remote sensors to validate the hypothesis. As the calculation directly depends on the NDVI values, one essential requirement is a previous, rigorous calibration and correction of the sensor reflectance to ensure a fair integration between the NDVIs from different dates. In this case, the radiometric calibration, together with the correction of the light conditions, guaranteed the stability of the images regardless of the illumination conditions and the sensor characteristics. Another remark that could be made is on the use of the red-edge band of the Micasense camera, profiting from its spectral capacity. However, this band was previously tested in comparison with crop parameters with the same dataset [44], and it was found that the indices based in the red-edge band did not correlate well with structural parameters such as LAI and biomass; however, it seemed better suited to the depiction of chemical parameters, as was also found in another similar research study [71].

The scale of observation is also a key factor to be accounted for. In this work, two scales of NDVI observations were used: point scale (approximately 1/8 m²) and plot scale (approximately 400 m²), both of which were enabled by the superhigh spatial resolution of the images. Other studies have compared regional vs. local results, profiting from the scalability of remote sensing images [58,65,66]. Predictions based on satellite platforms appear to improve with the increased radiometric purity of the pixels [59,65]. In this vein, UAS would seem a perfect alternative because the derived image resolution is much higher than that from satellite systems and because it excludes land cover mixed pixels and atmospheric effects. Our estimations using the UAS resolution were superior in terms of biomass and errors at the point scale, probably because the NDVI values that were strictly bound by the sampling area better fit the field estimations. In fact, when the NDVI was averaged onto a plot scale, the biomass variability between plots was slightly lost, particularly for the pea-based associations (Table 5). However, other authors [53] found that an integrated average NDVI determined using a window size larger than a 1 × 1 pixel improved the results. Considering that regional yield or biomass statistics usually provide an averaged value over large areas, it seems reasonable that the averaged values of the image may better fit those broad statistics. Further applications must then be applied at regional scales to appraise the method’s performance, particularly when using remote sensing sensors such as the Landsat series or Sentinel-2 or recent high-resolution systems such as Geosat-2. These optical sensors may allow long temporal NDVI series to be readily available, affording repeatability and continuous observations instead of labored drone campaigns.