Cover Crop Types Influence Biomass Estimation Using Unmanned Aerial Vehicle-Mounted Multispectral Sensors

1. Introduction

Cover crops in cropping systems are grown between cash crop seasons to incorporate into the soil or left as mulches before maturity. They provide several ecosystem services, including reduced soil erosion, improved soil health, carbon sequestration, weed suppression, and increased crop yield [1,2,3]. However, most of these benefits from cover crops depend on the amount of biomass they produce at termination [2,4,5,6].

The quantity of cover crop biomass produced determines the extent of the ecosystem services. For example, increased biomass of cover crops provides more soil cover, thus reducing soil erosion [7,8]. Weed suppression from cover crops increases with increasing biomass production [5,6]. Similarly, increased biomass of cover crops upon incorporation into the soil increases soil carbon sequestration [9,10]. Therefore, quantifying the aboveground biomass is essential, as it also estimates the belowground biomass and overall biomass production of cover crops [11,12].

The conventional approach to quantifying cover crop biomass involves hand sampling from a few spots and measuring the dry mass to estimate the biomass of the entire field. This destructive method is laborious, time-consuming, and subject to error. Therefore, the use of advanced non-destructive techniques is necessary for quantifying and mapping the aboveground biomass of large agricultural fields [13,14]. Moreover, using the shoot-to-root ratio of specific crops, we can estimate and model the overall biomass of different cover crops.

Several remote sensing techniques can be used to model the cover crop biomass. For example, multispectral images from Landsat-8 showed that vegetation indices based on red, green, blue, and near-infrared (NIR) reflectance can be used to estimate oat-radish biomass [15]. On the other hand, short-wave infrared-based vegetation indices from Sentinel-2 have been shown to have better estimation accuracy for various cover crop biomass types [16,17]. However, the low spatial resolution of satellite images makes it harder to estimate biomass precisely, especially for smaller-scale agricultural fields. On the other hand, high-resolution multispectral images from unmanned aerial vehicles (UAVs) can help address this issue in precision agricultural management.

One of the major advantages of using UAVs for crop phenotyping is the ability to utilize crop canopy height models (CHMs) and achieve high spatial resolutions [18]. Several studies have found promising results from CHMs on UAVs in estimating different cover crop biomass. However, their efficiency varies widely depending on the crop types and growth stages [19,20,21,22]. These conflicting findings suggest a need for improved crop species performance of CHMs and their utility in modeling biomass using vegetation indices.

The common vegetation indices used for biomass estimation have been unsuccessful due to index saturation at the later growth stages of the crops [14,23,24]. As the canopies cover the field, it becomes difficult to quantify the biomass using NIR band reflectance [18,25]. Moreover, the performance of these vegetation indices depends on the type of crop [22,26,27,28]. Therefore, crop-specific investigations are needed to understand the performance of different vegetation indices and CHMs in estimating biomass. Moreover, it is unknown how combining spectral signatures with CHMs can be used to model the biomass of different cover crops.

In this study, we used multispectral sensors mounted on a UAV to estimate the biomass of brassica, legume, grass, and three-species mix cover crops. These species have vastly different canopy structures. Therefore, identifying species-specific models is crucial for estimating their biomass and overall ecosystem services. Often, these species are mixed as a cover crop for multiple benefits, such as nitrogen fixation from legumes and deeper root systems from grasses [29,30]. Our goal was to use different band reflectance, vegetation indices, and CHMs to model the biomass of cover crops. Our specific objectives were (1) to identify the best model to estimate cover crop biomass at different growth stages from bitemporal multispectral images and (2) to study the influence of cover crop types on the model performance of biomass estimation through multispectral images.

2. Materials and Methods

2.1. Study Site

The study was conducted at the Texas A&M University Research Farm near College Station, Texas (30°33’6.69” N, 96°25’57” W, 68 m above mean sea level). The study site is located on the Brazos River floodplains and within the East Central Texas Plains ecoregion (Figure 1). The site’s climate is characterized as humid subtropical, with an average annual precipitation of 1060 mm and an average annual temperature of 20.6 °C [31]. The soil at the site is classified as Belk Clay (fine, mixed, active, thermic Entic Hapluderts) with a 0-1% slope and a silty clay or clay texture in the A and B master horizons [32].

The cover crops were established as part of an organic cotton field experiment initiated in 2020. The field study employed a randomized split-plot design, with tillage practices as main plots (conventional and minimum tillage) and cover crop management as subplots, replicated three times. However, the effect of the tillage practices was not considered in this study. Therefore, each cover crop type was replicated six times in the experimental design. Each experimental plot was 50 m × 8 m in size. Details of the field experiment can be found in Salehin, Rajan [30].

Cover crops planted in the cotton field were Austrian winter pea (Pisum sativum), which is a legume; purple-top turnip (Brassica rapa subsp. rapa); and oats (Avena sativa L.), both of which are non-legumes. A mix of all three cover crops was also compared with the monoculture cover crops. During the 2021-2022 season, the cover crops were planted on November 19, 2021, and terminated on April 6, 2022, after the final biomass sampling on April 1. The cover crop biomass samples were collected on March 1 and April 1, 2022, using two 0.5 × 0.5 m quadrants from each plot, a day after the flights. The two samples from each plot were merged to take the average for each plot. The biomass samples were oven-dried at 60 ºC for at least three days or until a constant mass was observed. The dry mass of samples from each plot was recorded for use in modeling with UAV images.

2.2. Image Acquisition

The images were collected using a five-spectral-band camera from MicaSense (RedEdge-3, MicaSense Inc., Seattle, WA). The five spectral bands, along with their central wavelengths and bandwidths, are presented in Table 1. The camera system utilized a 5.4 mm focal length and a 4.8 mm x 3.36 mm complementary metal oxide semiconductor (CMOS) sensor (Figure 2A). The camera was mounted on a UAV or quadcopter (Matrice 100, SZ DJI Technology Ltd., Shenzhen, China) (Figure 2B). The images were collected on February 28 and March 31, 2022, at the mid and late-growth stages of the cover crops.

The flights were conducted at a height of 30 m above ground with 75% front and side overlaps. The spatial resolution, or ground sampling distance, of the images was 2 cm, meaning that each pixel covered a 2×2 cm² area on the ground. The UgCS flight planning app (SPH Engineering, Riga, Latvia) was used to plan flight missions using an automated route. Four ground control points (GCP) were laid on four corners of the field in a zigzag pattern before the flight for accurate georeferencing for image processing (Figure 3). GPS coordinates and elevations of the GCPs were collected using real-time kinematics (RTK) system (Reach RS+, Emlid Ltd., Saint Petersburg, Russia) configured with Trimble network (Figure 2C).

2.3. Image Processing Approach

A total of 315 images were collected for each band on February 28, and 312 images were collected on March 31. Imagery for each date was processed using Pix4Dmapper (version 4.9.0, SA, Switzerland) in three major steps. In the first step, identifiable features were extracted as key points. Following this, the images were calibrated based on the camera model and geolocation of GCPs. In the second step, densified point clouds were produced using the tie points from the first step. The orthomosaic image was the product of the orthorectification process, which corrects lens distortion, camera tilt, perspective, and relief displacement [33]. The band reflectance values were adjusted based on calibration panel images captured at the start and end of each flight. The known reflectance values of the panel for each band are used to calibrate the camera reflectance. The final result was an orthomosaic dataset of single-band reflectance of blue, green, red, red edge, and NIR.

Individual plot shapefiles were created for each image acquisition date using ArcGIS Pro [34]. The overall statistics of individual band reflectance, including mean, were extracted using the plot shapefiles and the zonal statistics function in ArcGIS. The five single-band images were layered into a composite band for each date. The composite band was used to generate different vegetation indices using the band arithmetic tool in ArcGIS (Figure 4). The details of the calculated vegetation indices are provided in Table 2.

A rough triangular distribution of points was observed in the scatterplot when NIR reflectance was plotted against the corresponding red band reflectance (Figure 5A,C). This pattern is associated with the reflectance characteristics of vegetation and soil and the mixing of these characteristics within the image pixels. The diagrammatic triangular representation shows the features of the images in the distribution (Figure 5B,D). The bare soil areas within the images lie along the diagonal base of the triangle (Figure 5B,D). In contrast, pixels representing living vegetation are positioned in the upper-left region of the distribution. This distribution arises because living vegetation strongly absorbs light in the red spectral band while reflecting strongly in the NIR band [42]. However, the pixel distribution for soil and the full canopy does not always form straight lines, as shown in Figure 5B and Figure 5D. The intercept (a) and slope (b) of this bare soil line for each imagery are considered to calculate the perpendicular vegetation index (PVI).

The canopy height models (CHMs) were generated using Structure from Motion (SfM) point clouds created in Pix4Dmapper. Pix4D uses overlapped images to detect key points across multiple images. It then matched these key points to determine their relative position based on X, Y, and Z axes, which represent latitude, longitude, and elevation, respectively. Pix4D also outputs a digital surface model based on the densified point clouds. For further processing, the point cloud data was imported into Quick Terrain Modeler [43]. There, the original Z values (elevation) were then converted into the above-ground level (AGL) by subtracting the digital surface model from the original Z values using the AGL analyst tool. This AGL dataset was then used to create a gridded surface model representing the CHM. The February and March CHMs were generated at sampling distances of 11 and 10 cm (Figure 4B,D). Canopy height statistics were extracted from the CHMs using the shapefiles and the “Zonal Statistics as Table” function of ArcGIS. We considered the mean values of zonal statistics because of their high correlation with biomass. The detailed steps of image processing are shown in Figure 6.

2.4. Statistical Analyses

The mean values of the individual band reflectance, vegetation indices, and CHMs extracted from the plot shapefiles were used to model the biomass of cover crops. All statistical analyses were performed in R version 4.4.1 [44]. First, the correlation among the variables, including the cover crop biomass, band reflectance, vegetation indices, and CHMs, was investigated using the CORRPLOT package in R [45]. Later, simple linear regression models were fitted for each vegetation index and CHM to estimate biomass.

To determine the best-fit model, bidirectional stepwise multiple linear regression was performed separately for each date and cover crop type, based on P-values. For that, the biomass was modeled with all reflectance bands, nine vegetation indices, and the CHM. Then, both forward and backward stepwise regression were performed, using a P value of 0.15 to enter and remove variables from the model. The stepwise regression was performed using the OLSRR package in R [46]. The adjusted R² and root mean square error (RMSE) of these models were assessed to determine their performance. Due to the high correlation among some of the vegetation indices, we examined the variance inflation factor (VIF) for each model using the CAR package in R [47]. If the VIF was greater than 15 for a variable, one or more variables were removed from the model based on their correlation with biomass, thereby reducing the VIF of all variables without significantly compromising the model fit.

3. Results

For February images, all vegetation indices were highly correlated with the aboveground biomass (Table 3). Among the band reflectance, NIR had a high correlation coefficient with the biomass, followed by the red edge. On the other hand, the reflectance of blue, red, and green bands was negatively correlated with biomass. The modeled canopy height (CHM) did not correlate strongly with vegetation indices in February images (Table 3).

Similarly, in March images, vegetation indices were highly correlated with biomass. However, the correlation coefficients were not as high as in the February images. They ranged from 0.76 with NIR band reflectance to 0.91 for GNDVI and NDRE. Unlike in February images, the red edge band reflectance correlated negatively with the biomass in March images. However, the CHM had a better correlation with the biomass in March (Table 3).

When both image acquisition dates were considered, vegetation indices and band reflectance did not perform as well as when individual image acquisition dates were examined (Figure 7). NDVI, NDRE, PVI, and EVI showed the best correlation with biomass, followed by CIre, SRre, GNDVI, SR, CIg, and NIR reflectance. On the other hand, CHM has a correlation coefficient of 0.72 (Figure 7).

Simple linear regression models of vegetation indices and band reflectance with aboveground biomass were assessed to predict aboveground biomass. For February images, PVI, NDRE, SRre, EVI, and CIre had the highest R² values of 0.86 (Figure 8B–E,I). They were followed by SR, CIg, NDVI, and GNDVI with R² values of 0.85, 0.85, 0.84, and 0.83, respectively (Figure 8A,F–H).

Similarly, modeling cover crop biomass using vegetation indices extracted from the March images showed strong linear fits. For instance, NDRE and GNDVI had the best fit with an R² value of 0.82 (Figure 9A,B). Other vegetation indices, such as SRre, CIre, CIg, NDVI, EVI, SR, and PVI, have linear fits with R² values of 0.81, 0.81, 0.78, 0.76, 0.72, 0.69, and 0.73, respectively (Figure 9C–I).

On the other hand, when both February and March images were considered, NDVI had the best fit with an R² value of 0.79 (Figure 10A). This was followed by PVI with an R² value of 0.77 (Figure 10I). Other vegetation indices such as NDRE, EVI, SRre, CIre, SR, GNDVI, and CIg had linear fits with R² values of 0.76, 0.76, 0.74, 0.74, 0.72, 0.71, 0.71, and 0.73, respectively (Figure 10B–H).

The canopy height model from the February images didn’t have a strong relationship with the aboveground biomass. They had a linear fit with an R² value of 0.36 (Figure 11A). On the other hand, for images collected in March, the canopy height model had a quadratic fit with an R² value of 0.78 (Figure 11B). However, considering both image acquisition dates, the modeled canopy height had a quadratic fit with an R² value of 0.54 (Figure 11C).

Stepwise multiple linear regressions were performed to identify the best model fit for aboveground biomass. For this, a full model with five spectral band reflectance, nine vegetation indices, and the CHMs was fitted and analyzed for bidirectional stepwise regression based on P-values. The final model output for the February images included only PVI, with an R² value of 0.86 (Table 4, Figure 12A). In contrast, for the March images, NDRE, NDVI, and CHM were retained in the final model, yielding an R² value of 0.90. However, due to the high correlation between NDRE and NDVI, the variance inflation factor (VIF) was high for both NDRE (97.9) and NDVI (93.9). Therefore, based on the variable correlation with biomass, NDVI was excluded from the model. The final model, incorporating NDRE and CHM, had an R² value of 0.84 and an RMSE value of 408.3 (Table 4, Figure 12B). With the new model for March images, the VIF was reduced to 11.5 and 14.3 for NDRE and CHM, respectively.

When both image acquisition dates were considered, green band reflectance, CHM, SRre, and GNDVI were included in the final model, yielding an R² value of 0.86. However, due to the high correlation between the vegetation indices, the VIF was very high for SRre (63.2) and GNDVI (27.4). Therefore, GNDVI was removed from the model without compromising the model’s fit. The new model parameters had reduced VIF for Green (2.5), CHM (9.0), and SRre (11.9). The final model had an R² value of 0.85 and an RMSE of 345.8 (Table 4, Figure 12C).

Stepwise multiple linear regression analyses for each cover crop type showed that cover crop types interfere with the performance of models in biomass estimation. The full model considered the reflectance of all five bands, nine vegetation indices, and the canopy height model. The final model for oats had only EVI as the predictor, with an R² value of 0.86 and RMSE of 204 kg ha^-1 (Table 5, Figure 13A). The final model for AWP only had red edge band reflectance as the predictor, with an R² value of 0.71 and RMSE of 262 kg ha^-1 (Table 5, Figure 13B). However, turnips were best predicted by a combination of NIR reflectance, GNDVI, and CHM with an R² value of 0.95 and RMSE of 55 kg ha^-1 (Table 4, Figure 13C). Finally, the mixed species were predicted by the reflectance of NIR and blue bands with an R² value of 0.93 and RMSE of 179 kg ha^-1 (Table 5, Figure 13D).

4. Discussion

4.1. Temporal Differences in Biomass Estimation

In February images, all vegetation indices, including the NIR band reflectance, were highly correlated with cover crop biomass (Table 3). This indicates that the spectral signatures were better captured in the early to mid-growth stages. At this stage, the canopies are not fully closed. Therefore, the spectral signatures captured by the multispectral sensor can better correlate with the biomass at the early growth stages of crops [48,49]. During this period, vegetation indices such as NDRE, PVI, SR, EVI, CIre, SRre, and CIg can be used to model the biomass with similar performance (Figure 8).

However, later in the growth stages, as captured in March images, some vegetation indices fail to correlate with the biomass quantity. For example, the most popular index, NDVI, was not the best-performing index for modeling biomass (Figure 9F). Similarly, EVI, SR, and PVI also failed to model biomass confidently. This can be attributed to the AWP canopies covering the interrow spaces and capturing greater spectral signatures, but not necessarily correlating with biomass (Figure 9). Rather, using the red edge band in vegetation indices showed better performance in modeling biomass. Therefore, NDRE, SRre, and CIre performed better, followed by GNDVI (Figure 9A–D). This was influenced by strong negative correlations of biomass with green and red edge band reflectance in March images (Table 3). As the crop biomass increases in later growth stages, increased chlorophyll content can absorb more red-edge spectral energy [50,51]. Studies have reported the need for different vegetation indices depending on the growth stages of the crops [52,53].

When both image acquisition dates were considered, the overall performance of all vegetation indices and band reflectance was reduced. However, NDVI was still the best index, followed by PVI, NDRE, and EVI in modeling biomass (Figure 10A–C,I). This indicates that vegetation indices using NIR, red, red edge, and blue spectral signatures can effectively estimate biomass from multitemporal images. NDVI has been shown to correlate well with biomass across various crops, as it effectively captures the spectral signatures of plant chlorophyll content and differentiates between healthy and stressed vegetation [52,53,54].

However, NDVI can suffer from saturation issues in high biomass scenarios, which limits its effectiveness in densely vegetated areas [55,56]. We haven’t seen such a saturation effect in cover crop biomass estimation. The use of narrow-band vegetation indices, which incorporate the red-edge region, has also been emphasized. These indices can mitigate saturation effects and provide a more nuanced understanding of biomass dynamics [27,57,58,59]. As cover crops are usually terminated before maturity, NDVI and other vegetation indices can still serve as reliable predictors of biomass.

The modeled canopy height from images did not perform well in modeling cover crop biomass in the early growth stage (Figure 11A). In contrast, later in the growth stage, cover crop types contributed to a quadratic relationship between modeled canopy heights and biomass (Figure 11B). On the other hand, when both image acquisition dates were considered, the canopy height failed to model the biomass effectively. This was mainly due to the lack of association between the biomass and height of herbaceous AWP and broadleaf turnips. Kümmerer, Noack [20] found similar trends in the aboveground biomass of mixed species cover crops with CHMs. Only oat biomass had a good fit with their modeled canopy heights (Figure 11). Studies have found that grass species biomass can be better modeled by incorporating the canopy heights [59,60,61]. However, complex canopy structures can be challenging to model broadleaf or herbaceous crop biomass with canopy heights [61]. Moreover, the sensitive plant height range (2 to 7 cm) during early growth stages can influence model performance (Figure 11A).

As February images were well correlated with most vegetation indices, the stepwise regression found PVI to be the best fit in modeling cover crop biomass from early growth stages (Table 4, Figure 12A). PVI has the advantage of adjusting for bare soil reflectance. Therefore, it can capture vegetation traits more effectively than regular vegetation indices, especially during early growth stages when interrow bare soils are still exposed and captured in spectral signatures. However, later in the growth stage, spectral signatures get complicated, requiring multiple parameters to best model the biomass, as seen in March images (Table 4, Figure 12B). It required a combination of NDRE and CHM to model the biomass most effectively with an R² value of 0.84. Similarly, for bitemporal images, a combination of green band reflectance, CHM, and SRre was considered the best to model the biomass (Table 4, Figure 12C). This indicates that the use of multiple indices and CHM can improve the performance of cover crop biomass estimation in later growth stages and with multitemporal images [62]. This approach is particularly beneficial in heterogeneous environments where mixed pixels can complicate biomass assessments [63].

4.2. Influence of Cover Crop Types on Biomass Estimation

Cover crop types significantly impacted the performance of different indices and canopy height models in modeling biomass. For example, oats were best modeled using only EVI with a good fit (Table 5, Figure 13A). The narrow-band vegetation indices can correlate better with grass biomass, overcoming the saturation issues of NDVI [64]. Oats usually have significantly exposed soil in the interrow spaces compared to herbaceous crops. EVI incorporates blue band reflectance, which allows for adjustment of background soil reflectance. Therefore, it enhances the sensitivity to variations in biomass, particularly in high-density vegetation scenarios [65].

Similarly, AWP was best modeled by red edge band reflectance (Table 5, Figure 13B). Unlike the red bands used in NDVI, red edge band reflectance is sensitive to canopy structure in dense vegetation, such as AWP [66]. Due to the herbaceous canopy structure, regular vegetation indices focused on the NIR band were unable to accurately model the AWP biomass. Red edge band reflectance is typically lower than near-infrared (NIR) reflectance and is sensitive to high chlorophyll content and dense vegetation. Therefore, red edge reflectance can better model the dense AWP biomass.

However, turnip biomass was modeled using a combination of NIR, GNDVI, and CHM (Table 5, Figure 13C). Although individual indices and CHM did not perform well for estimating turnip biomass, this combination overcame the challenges and exhibited strong prediction performance with an R² value of 0.95. Corti, Cavalli [67] reported improved performance of herbaceous cover crop biomass modeling by incorporating UAV-derived CHMs with vegetation indices. However, the NIR reflectance and CHM had a negative relationship with the turnip biomass estimation in this model. This can be due to the spectral signatures and CHM being affected by the weed pressure in low-biomass turnip plots. The combination of NIR, GNDVI, and CHM can overcome this challenge and accurately model the aboveground biomass of the turnip.

A combination of NIR and blue band reflectance best modeled the biomass for mixed species (Table 5, Figure 13D). The strong NIR reflectance and absorbance of blue bands by mixed species can be effectively used to model biomass. The complex canopy structure of mixed species, including oats and AWP, makes it more challenging for vegetation indices to model biomass accurately. However, the use of band reflectance from multispectral sensors can be a valuable tool for modeling mixed-species biomass. Several studies have found the efficiency of vegetation indices in estimating legume-grass mix biomass with moderate accuracy (R² values ranging from 0.74 to 0.92) [68,69]. However, our study shows that using only the NIR and blue bands, we can accurately model oat, AWP, and turnip mix biomass with an R² value of 0.93.

5. Conclusions

This study represents the effectiveness of UAV-mounted multispectral sensors in estimating cover crop biomass, utilizing a combination of band reflectance, vegetation indices, and CHM. The findings indicate that vegetation indices, such as NDVI, NDRE, EVI, and PVI, performed well in biomass estimation, with varying efficacy depending on crop type and growth stage. In the early growth stages, spectral indices were strongly correlated (r = 0.91-0.93) with biomass, whereas in later growth stages, CHMs were required to integrate to improve model accuracy (R² = 0.84).

Additionally, the study highlights the need for different vegetation indices and band reflectance to model different cover crop biomass. Grass biomass, such as oats, can be best estimated by EVI (R² = 0.86), while legume species AWP can be estimated by red-edge reflectance (R² = 0.71). The low biomass and weed pressure in turnip plots necessitated a combination of NIR, GNDVI, and CHM, yielding an R² of 0.95. Finally, the complex canopy structure of mixed species can be best modeled by NIR and blue band reflectance (R² = 0.93).

A limitation of this study is the number of samples or temporal frequency needed to validate the model’s performance. For future investigations, additional experimental plots or sampling points should be included to test machine learning models and enhance their predictive performance. However, including multiple cover crop species allowed us to test the model for each species. The multiple linear regression models can estimate the biomass of oats, Austrian winter peas, turnips, and mixed species. Moreover, the known shoot-to-root ratios of each cover crop can be used to estimate the overall biomass. These findings also highlight the use of UAVs and multispectral sensors in precision agriculture. They can facilitate precise mapping of carbon and nitrogen capture, as well as other agroecosystem services, in large agricultural fields using cover crops.