Evaluation of Machine Learning Regression Techniques for Estimating Winter Wheat Biomass Using Biophysical, Biochemical, and UAV Multispectral Data

Abstract

Crop above-ground biomass (AGB) estimation is a critical practice in precision agriculture (PA) and is vital for monitoring crop health and predicting yields. Accurate AGB estimation allows farmers to take timely actions to maximize yields within a given growth season. The objective of this study is to use unmanned aerial vehicle (UAV) multispectral imagery, along with derived vegetation indices (VI), plant height, leaf area index (LAI), and plant nutrient content ratios, to predict the dry AGB (g/m²) of a winter wheat field in southwestern Ontario, Canada. This study assessed the effectiveness of Random Forest (RF) and Support Vector Regression (SVR) models in predicting dry ABG from 42 variables. The RF models consistently outperformed the SVR models, with the top-performing RF model utilizing 20 selected variables based on their contribution to increasing node purity in the decision trees. This model achieved an R² of 0.81 and a root mean square error (RMSE) of 149.95 g/m². Notably, the variables in the top-performing model included a combination of MicaSense bands, VIs, nutrient content levels, nutrient content ratios, and plant height. This model significantly outperformed all other RF and SVR models in this study that relied solely on UAV multispectral data or plant leaf nutrient content. The insights gained from this model can enhance the estimation and management of wheat AGB, leading to more effective crop yield predictions and management.

1. Introduction

With the increasing growth of the global population, the strong demand for food sources and food security has highlighted the need to enhance the development of efficient and sustainable agricultural practices. In the modern era, challenges such as rising global food demand, crop diseases and pest outbreaks, limited cultivated areas, and climate change are affecting the entire agriculture industry. Tan and Reynolds found that in southern Ontario, water supply and demand are the major challenges for the agricultural industry [1]. Notably, farmers in the province are less concerned about climate change compared to those in regions where extreme weather events are more prevalent [2]. The agriculture and agri-food sectors contribute approximately 7% to Canada’s gross domestic product (GDP), and one in every nine jobs in Canada was provided by this sector in 2022 [3]. Although climate change is not an immediate challenge for the Canadian agricultural industry, it is prudent to be informed early and prepare for counteractions while we still have time to respond to unforeseen climate variations.

Precision agriculture (PA) utilizes advanced technologies and data analysis techniques to maximize crop output while minimizing inputs. This approach involves assessing quantified spatial and in situ plant data to guide agricultural practices, such as the application of water, labor, and fuel, thereby reducing costs and preventing excessive waste, like pesticide and nutrient loss. PA integrates various spatial technologies, including geographic information systems (GIS), handheld ground-based data collection devices, and remote sensing through ground-based or aerial vehicles, to develop and implement efficient agricultural strategies [4].

Above-ground biomass (AGB) is a frequently used parameter to indicate crop growth status and the effects of agricultural management practices, making AGB estimation one of the main applications in PA [5,6]. In this study, we adopt a multivariate approach to estimate AGB using biophysical and biochemical parameters, utilizing in situ field data and high-resolution multispectral imagery collected by an unmanned aerial vehicle (UAV). Biophysical parameters included plant height and the leaf area index (LAI). In the context of using plant height as a predictor of AGB, the literature includes UAV-based height extraction methods that provide comprehensive coverage of the studied field, with multispectral cameras and LiDAR systems being common approaches [6,7,8]. Research has yielded varying degrees of success in identifying plant height as an important factor correlating with AGB. Furthermore, the LAI has been proven to be a significant parameter for monitoring crop growth and estimating AGB. Liu et al. found a strong linear relationship between the LAI and AGB, though this relationship weakens after the crops’ senescence [9]. To explore the potential of variables that are strong predictors of AGB in the early growth stages of winter wheat, our study attempts to address the limitations of these variables in later growth stages by incorporating biochemical parameters.

While research related to biomass estimation is abundant, few studies have utilized biochemical parameters, such as plant nutrient contents, as predictors. Common variables in biomass and yield estimation included plant height, the LAI, and specific vegetation indices (VIs), like the renormalized difference vegetation index (RDVI) and the modified hyperspectral variant of the normalized difference vegetation index (NDVI-like) [7,10,11]. These variables measure the physical structures of plants and are indicative of plant biomass and yield, and it is well established that plant biochemistry is intricately linked to plant structure, health, and condition, all of which are critical factors in applying PA strategies [12]. According to Marschner (2001) [13], macronutrients, micronutrients, and beneficial elements are essential classes of nutrients that promote plant health and growth through various mechanisms. For instance, nitrogen is a crucial macronutrient and is a major constituent of organic materials such as enzymes, chlorophyll, and compounds involved in oxidation-reduction reactions. The nitrogen content in plant tissue can indicate yield potential and overall crop health. Micronutrients, including iron, manganese, copper, and zinc, along with beneficial elements, like sodium, boron, and aluminum, play vital roles in plant growth. These micronutrients are essential for redox reactions and other physiological processes. For example, iron is necessary for protein synthesis and increases ribosome abundance in leaf cells. Manganese and copper act as activators for various enzymes, including those involved in detoxifying superoxide radicals and synthesizing lignin. Zinc is important for maintaining membrane integrity, protein synthesis, and the production of the phytohormone indole-3-acetic acid (IAA). Although beneficial elements are essential only for certain plant types, they stimulate growth and enhance physiological functions. Sodium, for instance, facilitates the movement of substrates between the mesophyll and the bundle sheath and can partially substitute for potassium’s role as an osmoticum. Boron contributes to cell wall stability by bridging polyuronides and promoting lignin synthesis. Understanding the significant roles of these nutrients underscores the potential of using plant nutrient contents as predictors of AGB. This approach could provide more comprehensive insights into crop health and productivity, thereby advancing the efficacy of PA practices.

The availability of plant nutrient data provides an opportunity to evaluate plant nutrient content ratios as predictors as well. Balanced nutrition is crucial for achieving high yields, and the overapplication of fertilizers can lead to reduced yields, soil and groundwater contamination, and harmful effects on human health and the environment [14,15]. While few studies have explored using plant nutrient content ratios as predictors of AGB, ratios, such as nitrogen to phosphorus (N:P), have been used in crop fertilization as indicators of nutrient limitations, particularly when either nitrogen or phosphorus is the limiting factor for plant growth [16]. The lack of research in this area, combined with the availability of relevant data, presents a promising opportunity to investigate the effectiveness of nutrient content ratios in estimating AGB.

Unmanned aerial vehicles (UAVs) are widely utilized in PA to capture timely, accurate, and cost-effective data on the earth’s surface [17]. Passive sensors, such as multispectral or hyperspectral cameras, RGB cameras, and active sensors, such as LiDAR, are typically mounted on UAVs to collect data for remote sensing applications in PA. These sensors are adopted because they do not require physical or destructive contact with plants to gather information. With the spectral data collected from remote sensing imagery, vegetation index (VI) calculations are made possible. VIs are mathematical transformations of spectral bands widely used in agricultural research to determine specific plant properties, such as the LAI, chlorophyll content, and nutrient levels [11,18,19]. Consequently, VIs are commonly adopted for crop growth and health monitoring, including biomass estimation, and research has demonstrated that VIs can be effective predictors of biomass [20,21]. For instance, vegetation indices that performed well in the study by Fu et al. were derived using the red absorption portion (550 nm–750 nm) of the spectrum [22]. On multispectral cameras, this typically includes the red band and red-edge bands. VIs that were proven by them as reliable predictors, such as the normalized difference vegetation index (NDVI) and the soil-adjusted vegetation index (SAVI), utilize spectral information from the red absorption portion. Based on these findings, it is imperative to further explore the biomass estimation capabilities of a diverse range of vegetation indices.

Although crop monitoring has traditionally relied on satellite imagery, UAV-based imagery offers significant advantages in terms of spatial and temporal resolution [23,24]. UAV systems can produce data with spatial resolutions of less than 10 cm compared to the meter-level spatial resolutions of satellite imagery. For example, the Landsat 8 Operational Land Imager, launched in 2013, has a spatial resolution varying between 15 and 30 m across its nine spectral bands and revisits the same geographical location every 16 days. Similarly, Sentinel-2, launched in 2015, features 13 multispectral bands with spatial resolutions of 10 m, 20 m, and 60 m, and a revisit time of 5 days with its constellation of twin satellites. Studies have indicated that UAV-based spectral data collected over smaller sampling areas explain more variation in wheat grain yield than the best-performing Sentinel-2 data. In comparison, Sentinel-2 data have yielded unsatisfactory results due to cloud coverage and lower temporal resolution [25]. This underscores the superior spatial and temporal resolution advantages that UAVs have over satellites for crop monitoring.

Winter wheat was selected for this study due to its prominence as one of the most widely cultivated crops in southern Ontario [26]. In recent years, machine learning regression methods, such as Random Forest (RF) and Support Vector Regression (SVR), have been extensively explored in biomass and yield estimation studies [27,28,29]. A significant advantage of these machine learning regression methods over linear regression is their applicability to a wide range of data, as they do not assume linear relationships. Given the diverse categories of variables involved, it is crucial to use a method suitable for capturing complex, non-linear relationships to ensure the validity of the results and reduce variability [30].

To make a well-informed estimation of AGB, it is essential to incorporate a wide range of data, including both biophysical and biochemical parameters. The objective of this study is to (i) investigate the relationships between AGB and factors, such as plant height, LAI, multispectral bands, VIs, and plant nutrient content levels and ratios; (ii) evaluate the effectiveness of RF and SVR models in estimating AGB; (iii) determine the optimal combinations of dates (growth stages) for AGB estimation in a winter wheat field located in southern Ontario; and (iv) identify the ranked importance and optimal combinations of variables for AGB estimation.

2. Materials and Methods

2.1. Study Area and Data Collection

The study site is in Southwest Middlesex County, Ontario, Canada, near the community of Melbourne, which is about 40 km southwest of the urban center of London, Ontario (Figure 1). Fieldwork was conducted in June of 2022, during which the average temperature was recorded at 18.8 °C and the relative humidity averaged 68.1%. The climate in the area is classified as warm summer humid continental climate (Dfb) according to the Köppen climate classification system. The area is predominantly agricultural croplands, and its major field crops include winter wheat, corn, and soybeans. Winter wheat was selected as the focus of this study. A winter wheat field covering 35.5 hectares (approximately 355,000 m²) in this region was designated as the specific area for investigation.

The cultivar in the studied field was soft red winter wheat, which was planted in October 2021. In the region of southwestern Ontario, winter wheat typically commences shooting in late April and is harvested from early to mid-July. Data acquisition was performed during the start of inflorescence emergence and heading stage to the ripening stage of the winter wheat. Studies have also pointed out that as AGB increases with the advancement of the growth stages, VIs’ correlation with AGB decreases after the flowering stage [28,31]. This decline is attributed to the maturation and yellowing of the plant’s leaves, underscoring the significance of incorporating more variety of data in assessing the efficacy of machine learning regression models in AGB estimation.

As outlined in

Table 1. Number of sample points and dates of data collection season.

Fieldwork Dates	Field Sample Point Groups	# of Sample Points	UAV Flight Dates	Phenology (BBCH Scale 1)
June 4	Winter Wheat Field W4 and W5	12 in W4, 16 in W5	June 8	Inflorescence emergence, heading (high 50 s to low 60 s)
June 10			June 10	Flowering, anthesis (60 s)
June 17			June 19	Development of fruit (70 s)
June 23			June 24	Ripening (low 80 s)

¹ Biologische Bundesanstalt, Bundessortenamt and CHemical industry.

, the field was revisited every six to seven days for ground sampling to align with the changing phenological stages of the winter wheat. Ground sampling was scheduled to align with UAV flights, whenever possible, to maintain consistency in data collection relative to the growth stages. However, optimal conditions for ground sampling and UAV flights were not always synchronized due to potential adverse weather conditions, such as strong winds or rapid weather changes. Given the field’s longest edge exceeding 850 m, two separate UAV flights were necessitated during each visit. Within the map depicted in Figure 2, 2 sets of sampling points were established: 16 sample points in a 4 × 4 grid on the northwest side and 12 sample points in a 4 × 3 grid on the southeast side. This placement aimed to maximize area coverage while minimizing labor intensity. The sampling points were positioned approximately 60 m apart, both vertically and horizontally, to ensure a representative distribution of data. A minimum distance of 50 m from all roads and houses was maintained to mitigate potential outliers and minimize disturbance to local residents. A GPS device was employed to facilitate precise revisits of these sample points during subsequent fieldwork sessions.

Fresh AGB samples were destructively harvested in a 20 × 20 cm grid at each sample point and transported to A&L Canada Laboratories for immediate fresh weight determination on the day of collection. Subsequently, these samples were dried in an oven at 60 °C for 48–72 hours. After drying, the biomass was weighed, and the top leaves of the plants were analyzed to determine nutrient content levels using the A&L PT2 plant test. This test provided nutrient content levels expressed in percentages and parts per million (ppm), as well as both actual and expected nutrient content ratios for each sample. The expected ratio serves as a target value for farmers, aimed at enhancing plant quality; it remains consistent across the field but varies according to the growth stage. In contrast, the actual ratio, derived from the actual nutrient content percentages, varied from sample to sample. In this study, both nutrient content levels and the actual ratio were utilized as biochemical parameters for predicting AGB. For the purposes of this paper, the actual ratio will simply be referred to as “ratio,” as the expected ratio is not utilized in this study.

Additionally, the LAI and crop heights were measured at each sample point as physiological parameters for the machine learning models. A LI-COR LAI-2200C equipped with a 180° view cap was utilized to measure a single LAI value at each sample point as the canopy of the wheat field densified approaching crop maturity. At each sample point, a recording sequence was employed, consisting of four readings above the canopy and eight readings below the canopy near the plant roots, distributed evenly with four in one row and four in the adjacent row. Scattering corrections were applied as needed during the above canopy recording procedures, contingent on ambient lighting and sky conditions. Furthermore, six individual plant height values were measured within a 1 m radius of each sample point using a meter stick, from which an average height value was calculated for each point. During the height measurements, the plants were left undisturbed to maintain their natural posture.

2.2. UAV Imagery

The UAV employed in this study was the Da-Jiang Innovations (DJI) Matrice 100, equipped with a MicaSense RedEdge narrowband multispectral camera (MicaSense Inc., Seattle, WA, USA), which collected spectral information across various bands (Figure 3). All flights were scheduled between 10 a.m. and 2 p.m. under cloud-free or near cloud-free conditions to minimize illumination variability across the field. Flights were postponed and rescheduled to the nearest possible date if the weather conditions were suboptimal, ensuring alignment with ground data collection and the phenological stages of plant growth. Additionally, flights were conducted under the lowest possible wind conditions to reduce challenges in image mosaicking due to plant movement. The flight plan was designed using the Pix4Dcapture app, which allows the pilot to adjust flight settings dynamically. At both sections (W4 and W5) of the study field, flights were conducted at altitudes of 50–60 m above ground level, often at the upper limits of the UAV’s manufacturer-recommended wind speeds. To preserve the data quality and flight efficiency, the UAV was set to fly at speeds between 3 to 4 m/s, depending on the windspeeds of the day. As the winter wheat matured and increased in height, the plants exhibited greater sway. Consequently, to ensure the accuracy of the resulting orthomosaics, all flights were performed with 85% front and side overlap** image capture. The flight paths were executed in a zigzag pattern, aligning with the orientation of the crop rows to enhance the accuracy of the orthomosaics. The outputs were weekly generated MicaSense band orthomosaics with a spatial resolution of 4 × 4 cm. Unfortunately, the flight data collected in the first week experienced a four-day delay relative to the sampling date due to adverse flight conditions. On June 4th, thin clouds scattered across the sky led to the initial assessment of the flight data as inaccurate, followed by three days of intermittent showers or wind speeds too high for safe UAV operation.

2.3. UAV Image Processing

Pix4Dmapper (version 4.8.0) was employed to process the multispectral images collected, generating one orthomosaic image per band. Prior to each flight, the MicaSense camera was subjected to a radiometric calibration to ensure the accuracy of the reflectance data. This calibration involved positioning the camera above a MicaSense Calibrated Reflectance Panel to capture white reference images for each band, taking into account sensor influences and the scene’s illumination conditions at the location. These white reference images, along with the reflectance values provided by the manufacturer’s white board, were utilized in Pix4Dmapper for image calibration, enabling each of the five MicaSense bands to produce corrected reflectance data of the field. The process of Structure from Motion (SfM), utilized by Pix4Dmapper, stitches together all the individual images captured by the camera. The high-overlap** image capture settings established for the flights facilitated this process, enhancing the accuracy of the results. The output comprised five orthomosaic images, each representing different reflectance values across the bands for both sections of the study area.

2.4. Vegetation Indices

The orthomosaic images were used to calculate vegetation indices (VIs) in QGIS. In order to minimize GPS error in the weekly visit to the field sample points, the VI values were averaged within a 1 m radius of the sample point. Details of the camera bands are listed in

Table 2. Spectral bands of the MicaSense multispectral camera.

Band	Name	Band Range (nm)	Center Wavelength (nm)	Bandwidth (nm)
1	Blue	465–485	475	20
2	Green	550–570	560	20
3	Red	663–673	668	10
4	Red Edge	712–722	717	10
5	NIR	820–860	840	40

.

A total of 13 Vis were calculated, as listed in

Table 3. Vegetation indices to be tested in this study.

VI 1	Formula 2	Authors
ARVI	$\frac{\mathrm{NIR}-[\mathrm{Red}-1\times \left(\mathrm{Blue}-\mathrm{Red}\right)]}{\mathrm{NIR}+[\mathrm{Red}-1\times \left(\mathrm{Blue}-\mathrm{Red}\right)]}$	Kaufman and Tanre [33]
Cl_RE	$\mathrm{NIR}÷\mathrm{RE}-1$	Gitelson et al. [34]
EVI	$\frac{2.5\times (\mathrm{NIR}-\mathrm{Red})}{\mathrm{NIR}+6\times \mathrm{Red}-7.5\times \mathrm{Blue}+1 }$	Huete et al. [35]
GCVI	$\mathrm{NIR}÷\mathrm{Green}-1$	Gitelson et al. [36]
ISR	$\mathrm{Red}÷\mathrm{NIR}$	Fernades et al. [37]
MCARI	$[\left(\mathrm{RE}-\mathrm{Red}\right)-0.2\times \left(\mathrm{RE}-\mathrm{Green}\right)]\times \mathrm{RE}÷\mathrm{Red}$	Daughtry et al. [38]
MSAVI	$[2\times \mathrm{NIR}+1-\sqrt{{\left(2\times \mathrm{NIR}+1\right)}^2-8\times \left(\mathrm{NIR}-\mathrm{Red}\right)}]÷2$	Qi et al. [39]
NDRE	$(\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{E})÷(\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{E})$	Gitelson and Merzyak [40]
NDVI	$(\mathrm{NIR}-\mathrm{Red})÷(\mathrm{NIR}+\mathrm{Red})$	Rouse et al. [41]
OSAVI	$[1.16\times (\mathrm{NIR}-\mathrm{Red}\left)\right]÷(\mathrm{NIR}+\mathrm{Red}+0.16)$	Rondeaux et al. [42]
RDVI	$(\mathrm{NIR}-\mathrm{Red})÷(\sqrt{\mathrm{NIR}+\mathrm{Red} })$	Roujean and Breon [43]
RVI	$\mathrm{NIR}÷\mathrm{Red}$	Jordan [44]
SAVI	$[1.5\times (\mathrm{NIR}-\mathrm{Red}\left)\right]÷(\mathrm{NIR}+\mathrm{Red}+0.5)$	Huete [45]

¹ ARVI, atmospherically resistant vegetation index; Cl_RE, chlorophyll index red edge; EVI, enhanced vegetation index; GCVI, green chlorophyll vegetation index; ISR, infrared simple ratio; MCARI, modified chlorophyll absorption in reflectance index; MSAVI, modified soil-adjusted vegetation index; NDRE, normalized difference red edge; NDVI, normalized difference vegetation index; OSAVI, optimized soil-adjusted vegetation index; RDVI, renormalized difference vegetation index; RVI, ratio vegetation index; SAVI, soil-adjusted vegetation index. ² Blue, blue reflectance; green, green reflectance; red, red reflectance; RE, red edge reflectance; NIR, near-infrared reflectance.

. Indices such as the NDVI and SAVI have been previously validated as reliable predictors of winter wheat biomass [22]. Additionally, several of the VIs make use of spectral information in the red edge and near-infrared wavelengths, which have been demonstrated to correlate strongly with crop growth, health, yield, and the LAI [11,32]. The chlorophyll index red edge (CL_RE) has been established as an effective VI for predicting crop nitrogen content, which serves as an indicator of plant vigor and productivity.

2.5. Biochemical Parameters

In this study, 14 nutrient content levels and 8 derived nutrient content ratios were analyzed. Existing research has identified nitrogen and phosphorus as essential for protein synthesis, enzyme activities, and chlorophyll formation in plants [46,47]. Additionally, potassium is crucial in mitigating stress from drought, cold temperatures, salinity, and biotic factors, such as diseases and pests [48]. For example, sufficient potassium levels can enhance photosynthetic efficiency, improve water usage, and stabilize plant metabolism under drought conditions. Additionally, nutrient content ratios, such as nitrogen to sulfur (N:S) in plant leaves, are significant indicators of crop health and nutrient deficiency [49,50]. This framework provided the basis for testing both individual nutrient content levels and ratios. The 14 nutrients tested included nitrogen (N), phosphorus (P), potassium (K), magnesium (Mg), calcium (Ca), sodium (Na), sulfur (S), boron (B), zinc (Zn), manganese (Mn), iron (Fe), copper (Cu), aluminum (Al), and nitrate-N. The 8 nutrient content ratios evaluated were N:S, N:K, P:S, P:Zn, K:Mg, K:Mn, Ca:B, and Fe:Mn.

2.6. Machine Learning Regression Modeling

In the context of machine learning, regression models are used to predict continuous outcomes based on input variables. Two prominent techniques within this domain are Random Forest (RF) Regression and Support Vector Regression (SVR), both of which offer robust solutions to complex regression problems.

RF is an ensemble learning method that operates by constructing multiple decision trees during the calibration phase and outputting the mean prediction of the individual trees. This method capitalizes on the power of multiple decision trees to reduce overfitting, which is common in models relying on a single decision tree. Each tree in the forest is built from a random sample of the calibration data, and at each node, a subset of features is randomly chosen to decide the split. This randomness helps in making the model more resilient to noise in the dataset. Moreover, Random Forest can handle large datasets with higher dimensionality and can estimate which variables are important in the underlying relationships being modeled.

SVR, on the other hand, extends the concepts of Support Vector Machines (SVMs) from classification to regression. Unlike traditional methods that minimize the error between predicted and actual values, SVR attempts to fit the error within a certain threshold. It involves the creation of a hyperplane in a multidimensional space where the distance between the data points and the hyperplane is minimized, ensuring that errors do not exceed a defined threshold. This makes SVR particularly useful in cases where a margin of tolerance is specified in the predictions. SVR is highly effective in handling non-linear relationships through the use of kernel functions, which map input data into higher-dimensional spaces [51].

Both RF and SVR provide distinct advantages depending on the nature of the data and the specific requirements of the regression task [29]. RF is generally preferred for problems with high-dimensional spaces and large datasets, offering interpretations in terms of feature importance. SVR is advantageous when dealing with datasets where the prediction needs to stay within a certain range and is effective in capturing complex relationships through its kernel trick. When employed thoughtfully, both methods can yield highly accurate predictive models in a wide range of scientific and industrial applications.

Figure 4 displays the workflow of the methodology. The modeling was written in R programming language using R Studio by utilizing packages such as “randomForest” and “e1071” for RF and SVR, respectively. In both models, the independent variables were the VIs, MicaSense bands, plant physiological parameters, plant nutrient levels, and plant nutrient content ratios. Data collected over the four weeks were randomly divided into a 70% calibration set and a 30% validation set. The strength of the prediction model was assessed using the coefficient of determination (R²) and root mean square error (RMSE). To ensure the model’s strength and generalizability, we validated the results by creating random splits of the calibration and validation sets 100 times. The R² and RMSE values reported are the average values obtained from these splits. The equations for both metrics are as follows:

$${\mathrm{R}}^2=1-\frac{{\sum \left({\mathrm{y}}_{\mathrm{i}}-{\widehat{\mathrm{y}}}_{\mathrm{i}}\right)}^2}{{\sum \left({\mathrm{y}}_{\mathrm{i}}-{\stackrel{\mathrm{-}}{\mathrm{y}}}_{\mathrm{i}}\right)}^2}$$

(1)

where ${\mathrm{y}}_{\mathrm{i}}$ is the observed value, ${\widehat{\mathrm{y}}}_{\mathrm{i}}$ is the predicted value, and ${\stackrel{\mathrm{-}}{\mathrm{y}}}_{\mathrm{i}}$ is the mean of the observed values, and

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

(2)

where ${\widehat{\mathrm{y}}}_{\mathrm{i}}$ represents the predicted AGB (g/m²), ${\mathrm{y}}_{\mathrm{i}}$ denotes the observed AGB (g/m²), n is the total number of observations, and i serves as the summation index, incrementing by one.

3. Results

3.1. Biomass Data

AGB was destructively collected at each sample point. The dry weight of the sampled biomass progressively increased across different growth stages, as illustrated in Figure 5. Initially, AGB exhibited a modest increase during the first two weeks of fieldwork, followed by a significant acceleration in growth thereafter.

3.2. Regression Models with All Variables

A total of 42 variables were utilized as predictors for AGB, including plant height, LAI, MicaSense bands, VIs, and levels and ratios of plant nutrient content. The datasets were categorized into single-date and multi-date groups to assess the temporal impact on the models and to identify the most effective date or combination of dates for estimating AGB. These variables were incorporated into the calibration and validation of the RF and SVR models, as detailed in

Table 4. Calibration and validation statistics: analysis by date and modeling approach (RF and SVR) using 42 variables, including plant height, the LAI, MicaSense bands, vegetation indices, and plant nutrient content levels and ratios ¹. n is the number of data entries.

Date	Model	(n)	Calibration		Validation
			R2	RMSE (g/m2)	R2	p-Value	RMSE (g/m2)
June 4	RF	28	0.95	41.90	0.21	NS	137.37
	SVR	28	0.93	40.11	0.47	<0.05	132.17
June 10	RF	28	0.96	54.64	−0.13	NS	86.19
	SVR	28	0.85	58.74	−0.14	NS	99.63
June 17	RF	28	0.93	52.62	−0.02	NS	162.74
	SVR	28	0.76	76.20	−0.14	NS	132.08
June 23	RF	28	0.95	80.82	−0.14	NS	245.32
	SVR	28	0.70	113.80	−0.14	NS	238.80
June 4, 10	RF	56	0.95	44.03	−0.06	NS	134.58
	SVR	56	0.75	63.90	0.09	NS	130.63
June 4, 17	RF	56	0.94	60.62	0.57	<0.001	138.11
	SVR	56	0.85	77.20	0.47	0.001	155.47
June 4, 23	RF	56	0.97	75.91	0.68	<0.001	237.24
	SVR	56	0.96	83.87	0.59	<0.001	257.60
June 10, 17	RF	56	0.92	59.01	0.40	<0.01	123.23
	SVR	56	0.84	71.79	0.46	0.001	120.43
June 10, 23	RF	56	0.97	71.46	0.66	<0.001	212.96
	SVR	56	0.96	83.51	0.68	<0.001	201.40
June 17, 23	RF	56	0.94	86.27	0.23	<0.05	252.59
	SVR	56	0.90	103.30	0.22	<0.05	249.10
June 4, 10, 17	RF	84	0.94	58.61	0.41	<0.001	131.64
	SVR	84	0.84	73.80	0.38	<0.001	130.17
June 4, 10, 23	RF	84	0.96	82.96	0.67	<0.001	207.08
	SVR	84	0.94	98.08	0.71	<0.001	184.32
June 4, 17, 23	RF	84	0.94	91.83	0.76	<0.001	177.09
	SVR	84	0.90	112.88	0.72	<0.001	187.79
June 10, 17, 23	RF	84	0.94	89.50	0.72	<0.001	177.91
	SVR	84	0.90	108.79	0.77	<0.001	156.61
June 4, 10, 17, 23	RF	112	0.93	90.98	0.80	<0.001	152.71
	SVR	112	0.89	113.81	0.77	<0.001	165.71

¹ All calibration models are significant at p-value < 0.001.

. Overall, the RF models exhibited slightly superior performance compared to the SVR models, with multi-date RF models outperforming those based on single dates. The best-performing RF model, which utilized all variables across all four dates, achieved an R² of 0.93 and an RMSE of 90.98 g/m² in its calibration set, and an R² of 0.80 with an RMSE of 152.71 g/m² in its validation set. RF models that incorporated data from three dates also demonstrated high performance. Similarly, SVR models showed improved performance with multi-date data compared to single-date models. The optimal SVR model, employing all variables and data from June 10, 17, and 23, yielded an R² of 0.90 and an RMSE of 108.79 g/m² in its calibration set, and an R² of 0.77 with an RMSE of 156.61 g/m² in its validation set. Although this model was only marginally superior to its counterpart, which utilized data from all four dates, it featured a lower RMSE. It is noteworthy that almost all models based on single-date data were not significant, a result anticipated due to the high number of variables relative to the modest dataset size of 28 entries.

It is crucial to note that the overall best-performing model is not necessarily the one with the highest R² value in either the calibration or validation sets. For example, the RF model using data from June 10 and 23 demonstrated a high R² of 0.97 and a low RMSE of 75.91 g/m² in the calibration set. However, the same model exhibited significantly weaker performance in the validation set, with an R² of 0.66 and a high RMSE of 212.96 g/m². This discrepancy suggests potential overfitting, indicating that while the model predicts the calibration data exceptionally well, it does not generalize effectively to new, unseen data.

3.3. Variable Importance Plot

RF modeling, which involves the use of numerous decision trees, was employed to generate a variable importance plot in R Studio using the “varImpPlot()” function. The plot displays increasing node purity (IncNodePurity) on the x-axis, which indicates the importance of each explanatory variable in predicting dry AGB on the decision trees. A higher IncNodePurity value signifies that the variable is more critical as a predictor. This method was employed to visualize the variable rankings in both the RF and SVR models, aiding in the identification of key predictors in the models.

The RF model incorporating all 42 variables demonstrated optimal performance when applied to the full four-date dataset. Analysis of the variable importance plot revealed that the NDVI was the most critical predictor, as depicted in Figure 6. Among the top ten most influential variables, the composition included five of the thirteen VIs utilized, two out of fourteen nutrient content levels, two of the eight nutrient content ratios, and plant height. Notably, the NDVI, ISR, ARVI, and RVI exhibited significantly higher IncNodePurity values compared to the remaining variables. These VIs are commonly associated with vegetation monitoring in agriculture and biomass estimation. Macronutrients such as N and K were also ranked in the top 10. N and K are crucial macronutrients that regulate enzymes and the synthesis of organic compounds. Additionally, K plays a vital role in cell growth and the regulation of photosynthesis, both of which are responsible for plant development [13,48].

The SVR model that incorporated all 42 variables yielded the best results using data collected on June 10, 17, and 23. A variable importance plot generated from these three dates identified the ISR as the most crucial predictor of AGB, with a slightly higher ranking than the NDVI, as shown in Figure 7. The top ten most important variables included five of the thirteen VIs utilized, three of the fourteen nutrient content levels, and two of the five MicaSense bands. Indices such as the NDVI, ISR, ARVI, and RVI displayed significantly higher IncNodePurity values than the other variables. Consistent with the results from the four-date analysis, the primary predictors remained the ISR, NDVI, RVI, and ARVI, albeit in a different order. In this three-date plot, MicaSense red and green bands, along with P, saw an increase in their rankings, moving into the top ten, a shift from their positions in the four-date plot. P is an essential macronutrient similar to N and K. However, its relevance had only increased in the four-date plot, while the importance of N and K had decreased. P is responsible for various biochemical reactions within the plant, including nitrogen fixation and the synthesis of nucleic acids and phospholipids, making it essential for the genetic and structural components of plant cells [13]. Although KMg_ACT was the highest-ranked nutrient content ratio in both plots, none of the nutrient content ratios were among the top ten variables in the three-date analysis.

3.4. Regression Models with Selected Variables

Variable selection is essential for reducing redundancy and complexity in regression models with a larger variety of variables. It enhances model performance and interpretability by focusing on the most relevant predictors, thereby reducing overfitting and computational complexity. This process also improves the model’s ability to generalize well to new data, ensuring more robust and accurate predictions. As indicated by the variable importance plots (Figure 6 and Figure 7), although the best-performing RF and SVR models utilized all 42 variables, the importance of the explanatory variables in estimating AGB varied considerably. For the RF models, based on the variable importance plot in Figure 6, we selected the top four, seven, ten, fourteen, twenty, and twenty-nine variables, determined by significant reductions in the increase in node purity, which helped establish a ranking threshold. A similar selection process was applied to the SVR models, grou** variables into rankings of the top five, seven, ten, fourteen, twenty, and twenty-eight based on the variable importance plot in Figure 7.

Additionally, two distinct groups of variables were evaluated in both RF and SVR models. Research has substantiated that UAV multispectral data alone can effectively predict AGB, leading to the testing of a class consisting solely of MicaSense bands and VIs [52]. Moreover, plant nutrient content, though seldom used as a predictor for AGB, displayed high importance in some explanatory variables. Consequently, nutrient content levels and ratios were also explored as a separate class for testing.

Table 5. Statistics of the RF models for above-ground biomass estimation with all dates (June 4, 10, 17, 23) and different combinations of variables (n = 112) ¹.

Variables	Number of Variables	Calibration		Validation
		R2	RMSE (g/m2)	R2	RMSE (g/m2)
All VIs + 5 MicaSense bands	18	0.91	102.19	0.73	175.63
All plant nutrient content + ratios	22	0.93	100.79	0.68	196.54
Top 4: NDVI, ISR, ARVI, RVI	4	0.87	119.81	0.59	213.49
Top 7: top 4 + height, N, KMg_ACT	7	0.91	100.90	0.76	167.32
Top 10: top 7 + NS_ACT, K, OSAVI	10	0.93	92.36	0.79	156.00
Top 14: Top 10 + red, NK_ACT, GCVI, NDRE	14	0.93	89.53	0.78	160.04
Top 20: Top 14 + green, Cl_RE, P, Fe, RE, Mg	20	0.93	89.19	0.81	149.95
Top 29: top 20 + NO3N, Ca, PZn_ACT, Al, MSAVI, MCARI, RDVI, CaB_ACT, Cu	29	0.94	89.41	0.81	151.52

¹ All models are significant at p-value < 0.001.

outlines the statistics for RF model calibration and validation sets using data from various dates and combinations of variables. The most effective date combination for RF utilized data from all four dates. For the calibration sets, the R² remained high across most groups, consistently above 0.9, except for the group containing only the top four variables. The RMSE values ranged from 89.19 to 119.81 g/m². Notably, the RMSE values were significantly elevated for groups comprising solely multispectral data and VIs, as well as those limited to plant nutrient content levels and ratios, with both exceeding 100 g/m². A decreasing trend in the RMSE was observed as more variables were included in the calibration sets, continuing up to the top 20 variables. Although the calibration sets exhibited similar performance, the validation set using only multispectral data and VIs outperformed the set that included only plant nutrient content levels and ratios. Nevertheless, neither was the best performing among the RF models. The validation set with the top seven variables demonstrated higher performance, which continued to improve up to the model incorporating the top twenty variables. The validation sets exhibited R² values between 0.59 and 0.81, with RMSE values spanning from 149.95 to 213.49 g/m².

The overall best-performing RF model employed a combination of the top 20 variables, achieving an R² of 0.93 and an RMSE of 89.19 g/m² in the calibration set and an R² of 0.81 and an RMSE of 149.95 g/m² in the validation set. Analyses of the models with the top 29 variables and all 42 variables (as detailed in Table 4) indicated a decline in model performance with the addition of more variables. The R² values plateaued while the RMSE increased, suggesting that eliminating lower-ranked variables from the variable importance plot can enhance the performance of the RF models.

Table 6. Statistics of the SVR models for above-ground biomass estimation with the three dates (June 10, 17, 23) and different combinations of variables (n = 112) ¹.

Variables	Number of Variables	Calibration		Validation
		R2	RMSE (g/m2)	R2	RMSE (g/m2)
All VIs + 5 MicaSense bands	18	0.81	145.51	0.62	190.51
All nutrient content + ratios	22	0.85	136.75	0.69	206.40
Top 5: ISR, NDVI, RVI, ARVI, K	5	0.72	178.05	0.66	187.07
Top 7: Top 5 + P, Red	7	0.76	162.81	0.62	198.99
Top 10: top 7 + GCVI, Na, green	10	0.81	147.21	0.66	184.68
Top 14: top 10 + RE, Cl_RE, blue, NDRE	14	0.81	145.18	0.69	179.08
Top 20: top 14 + OSAVI, B, KMg_ACT, height, NO3N, KMn_ACT	20	0.86	128.30	0.73	165.47
Top 28: Top 20 + CaB_ACT, NIR, SAVI, MSAVI, NS_ACT, N, PZn_ACT, Fe	28	0.88	119.65	0.77	154.36

¹ All models are significant at p-value < 0.001.

presents the statistics for SVR models using data from June 10, 17, and 23 across various variable combinations. The calibration sets of the SVR models displayed a range of performance with R² values between 0.72 and 0.88 and RMSE values from 119.65 to 178.05 g/m². The highest performance was observed with the top twenty-eight variables and the lowest with the top five variables. There was a linear increase in R² and a corresponding decrease in the RMSE as the number of variables in the calibration sets increased. Unlike the calibration sets, the validation sets did not exhibit consistent trends, with R² values ranging from 0.62 to 0.77 and RMSE values between 154.36 and 206.40 g/m². In contrast to the RF models, the SVR model utilizing only multispectral data and VIs performed worse than the model focusing solely on plant nutrient content levels and ratios. However, neither of these models achieved the highest performance.

The best overall performing SVR model utilized the top 28 variables, achieving an R² of 0.88 and an RMSE of 119.65 g/m² in the calibration set and an R² of 0.77 and an RMSE of 154.36 g/m² in the validation set. Comparing the performance of the model with the top 28 variables to that using all 42 variables (as detailed in Table 4), the former exhibited slightly better generalization capabilities, as indicated by a lower RMSE in the validation set. Therefore, the SVR model with 28 variables was considered superior due to its slightly better balance between training accuracy and validation error. Nonetheless, the differences in validation performance were minimal, suggesting that both models are relatively comparable in their ability to generalize.

4. Discussion

In this study, RF and SVM regression methods were used to predict the AGB of winter wheat, utilizing UAV multispectral MicaSense bands, associated VIs, plant biophysical parameters (plant height and the LAI), and plant biochemical parameters (nutrient content levels and ratios). During the first two weeks of sampling, the variation in AGB was low; however, it began to increase rapidly starting from the third sampling date, June 17. This rapid change aligns with the growth stages of winter wheat. On June 4 and 10, the plants were in their late heading and early flowering stages, respectively—a transition period marked by a slowdown in height increase due to shifts in developmental priorities and physiological changes. Initially, plant height and leaf area are major contributors to AGB as the stem elongates and leaves enlarge. By the heading stage, most stem elongation is complete, with culms extended and plant height largely established. What follows is primarily the emergence of the inflorescence from the flag leaf’s sheath, which does not significantly contribute to further height increase. As the plant transitions to the reproductive phase, its focus shifts from vegetative to reproductive growth, including the formation and maturation of the inflorescence. On June 17, the winter wheat entered the fruit development stage, channeling photosynthates from leaves and stems into the develo** grains, which accumulate starch, proteins, and other nutrients, significantly increasing their weight and overall biomass. Finally, during the ripening stage on June 23, the grains transform from a watery, milky substance into a hard, dry state—a process marked by the continuous accumulation of dry matter, primarily starch, enhancing total AGB.

The RF and SVR models were initially calibrated using all 42 variables across single- and multi-date datasets. In the validation sets, models based on single dates generally performed poorly and were statistically non-significant. Our finding agrees with Atkinson Amorim et al.’s work that models using multiple date combinations performed better, especially if the combinations involved the last two dates [27]. The best performance was observed in the RF model that utilized data from all four dates. According to the corresponding variable importance plot, the NDVI emerged as the most crucial variable, with other top-ranked variables primarily consisting of multispectral data. This finding aligns with Fu et al.’s findings that the NDVI and its narrowband-modified variations were effective predictors of biomass using partial least squares regression. Moreover, the top-ranked VIs in AGB estimation mainly comprised the NIR band, further supporting Fu et al.‘s identification of NIR as a sensitive band region for AGB [22]. In PA, employing UAVs equipped with spectral cameras is a common and cost-effective method to capture multispectral data and estimate AGB [53,54,55]. Although previous studies have proposed new methods for estimating crop AGB using multispectral data, our study demonstrates that existing machine learning methods can also produce comparable results when increasing the variety of predictors, such as plant nutrient content levels and ratios. Furthermore, senescence might affect the use of multispectral data for crop biomass monitoring at post-flowering stages, potentially reducing model performance compared to pre-flowering stages [56]. Our study did not encounter this issue, likely due to the inclusion of additional variables beyond multispectral data and derived VIs.

We successfully optimized the machine learning models by selecting the top-ranked variables from the variable importance plots. This approach is consistent with findings from other research, which has proven that machine learning is an effective method for predicting biomass in crops, such as wheat and oats [27,28,56]. While we were using a very similar camera setup installed on the UAV as Sharma et al., our findings proved that including biophysical and biochemical parameters in the analysis can significantly increase the RF and SVR models’ accuracy [56]. Similar to the results of Lu et al., the best-performing RF model proved more accurate than the top SVR model [57]. We tested a total of 42 variables. Although the SVR model performed well, nearly matching the RF model, the latter was better in terms of performance and ease of use. In the RF model, the top twenty variables were selected, which included eight of the thirteen VIs, three of the five MicaSense bands, five of the fourteen nutrient content levels, three of the eight nutrient content ratios, and plant height as the model’s performance started drop** with more variables being added. In comparison, Wang et al. reported higher RMSE values in their post-flowering stage analysis than ours in their linear regression, partial least squares regression (PLSR), and RF models [28]. Our model performance was comparable with theirs, and the lower RMSE values in our models could be advantageous for making timely adjustments in fertilizer and water application recommendations. This is especially crucial in unstable climate conditions, where late-stage growth adjustments are necessary. The RF method, which utilizes multiple decision trees, seems particularly well-suited to handling a high number of variables and avoiding overfitting. As reported by Wang et al., RF outperformed linear regression and PLSR by having higher stability in predicting AGB during their study period. Our findings aligned with theirs, as RF proved to be the more stable model when handling a variety of data compared to SVR. Conversely, the SVR model required user hyper-tuning and used a kernel trick function to separate data into groups, relying on the radial distance between points to provide meaningful insights for the model. We believe that our research findings can be applied to North American wheat fields located at similar latitudes to southern Ontario, Canada. However, the model may perform differently under varying agricultural practices and environmental conditions. Future applications of these findings should take these factors into account.

5. Conclusions

5.1. Contributions of Utilizing Multiple Categories of Variables in AGB Estimation

This study tested the effectiveness of multispectral data and biophysical and biochemical parameters in predicting the AGB of winter wheat using machine learning methods. Variables tested include UAV MicaSense bands, the derived VIs, plant height, the LAI, and plant nutrient content levels and ratios. The best result was obtained from the RF model with an R² of 0.81 and an RMSE of 149.95 g/m² using the top 20 variables, with a close to even split between spectral variables and nutrient content variables.

The inclusion of plant nutrient content levels and ratios as predictors in this study represents an advancement in the field of biomass estimation. Traditionally, these variables are not commonly utilized. The utilization of a lower-cost UAV multispectral camera setup, combined with biophysical and biochemical parameters, particularly at the later growth stages of winter wheat post-flowering, demonstrates a cost-effective method to predict AGB when urgent changes in late growth stages are needed to counteract unpredicted weather events, such as forest fire smoke and haze.

5.2. Limitations and Future Work

Though machine learning algorithms are capable of analyzing variables across categories, it is important to recognize the empirical nature of the machine learning models used. These models, by design, rely on existing datasets for validation and can only approximate the true AGB, which is only verifiable at harvest. This intrinsic limitation highlights the potential discrepancies between predicted and actual outcomes. Such limitations underscore the necessity for ongoing calibration and testing of these models under varied agricultural conditions and across different crop cycles to ensure their reliability and accuracy. The applications of these models can also be limited to dataset access, which is a common limitation in PA research because in situ measurements often are required, and that comes with the associated intensive labor, costs, and conditions. Additionally, it is important to recognize that the variables tested are not the only associated factors that affect AGB. Variables such as weather conditions, soil properties, field topography, moisture supply, and more need to be considered to define the condition of the plants.

Due to constraints in data and time availability, this research could not be conducted earlier in the growth stages of the winter wheat. Future studies should consider extending the time span to investigate the models’ effectiveness more comprehensively. Additionally, further exploration into the use of UAV-based spectral data for biomass estimation is suggested. This exploration should particularly focus on wavelengths or VIs that strongly correlate with plant nutrient levels and ratios. More importantly, it should emphasize hyperspectral bands, which have proven to be highly accurate in monitoring crop growth and estimating yield [58]. Furthermore, high spatial and temporal resolution satellite imagery can serve as a viable alternative to UAV imagery, eliminating the need for UAVs. Examples such as Planetscope and VENμS both have frequent revisit periods and high spatial resolution for a local, field-scale study. We demonstrated the effectiveness of using plant nutrient content levels and ratios as parameters to estimate AGB in this study. Therefore, further research into non-destructive methods using remote sensing techniques to obtain these data is recommended for future biomass estimation studies. Since plant height has also been proven to be a reliable predictor in this study, integrating Real-Time Kinematic (RTK) UAVs or LiDAR-equipped UAVs could enhance the precision and quantity of height data collection across the entire field. Combining these technologies with biomass estimation models could lead to the development of highly accurate AGB estimation maps, providing a more detailed understanding of biomass and crop yield potential.