Combining multiple UAV-Based indicators for wheat yield estimation, a case study from Germany

ABSTRACT

Unmanned aircraft vehicles (UAV) are widely used for yield estimations in agricultural production. Many significant improvements have been made towards the usage of hyperspectral and thermal sensors. The practical application of these new techniques meanwhile has been limited by the cost of data collection and the complexities of data processing. The objective of this paper is to evaluate the effectiveness of wheat yield estimations based on integrating vegetation indices (VI), solar radiation and crop height (CH), all of which are characterized by lower cost of data collection and processing. The VIs, solar radiation and CH were calculated based on UAV-based multispectral images obtained from two separate plots in Southern Germany and validated with data from a third plot. We compare the individual and joint predictive performance of different VIs, CH, and solar radiation by contrasting the estimated yield with actual yield based on multiple linear regression and quantile regression. The best predictive power was found for a combined estimation with CH, solar radiation and a Normalized Difference Red-edge Index (R² = 0.75,  RMSE = 0.53). This combined estimation resulted in a 15–20% improvement in the prediction of wheat yield accuracy as compared with utilizing any of the indices separately.

KEYWORDS:

• Drone sensors
• crop height
• multiple linear regression
• quantile regression
• multiple indicators
• crop surface model

Introduction

Wheat has become one of the most expansively cultivated food crops in the world (Tao et al., 2020). Because of its high economic and nutritional value, monitoring wheat crop growth and yield parameters during the season is essential to decision-makers (Marino & Alvino, 2020). The prompt and accurate estimation of wheat yield is a vital component of precision farming, and plays a definitive role in shaping food security and price regulation, therefore enhancing the stability of the agrarian economy (Tao et al., 2020). However, traditional crop yield estimations are reliant on ground-based data collection, which requires large time and financial resources (Reynolds et al., 2000) and therefore is unsuitable for upscaling (Huang et al., 2016; Liu et al., 2017). In recent years, remote sensing via satellite or drone technology has increased the spectrum of options (Viljanen et al., 2018). This allows for real-time yield estimations and crop health monitoring. Unmanned Aircraft Vehicles (UAV, or drones) use photogrammetric imaging and are configurated from motion techniques and tight image combination, providing 3D information at high flexibility and resolution.

The suitability of UAV-based information for high-precision yield estimations has been established by a growing body of literature: Recently, Tao et al. (2020) analyzed crop yield estimation based on wheat plant height information and vegetation indices obtained by a drone-mounted hyperspectral camera. Their model incorporated plant height parameters and vegetation indices obtained from UAV-based hyperspectral images. Similarly, Possoch et al. (2016) used multi-temporal crop surface models extracted from UAV hyperspectral data for biomass monitoring of barley. Their results reveal that plant height information combined with vegetation indices has the highest correlation with dry biomass.

The selection of actual prediction parameters, meanwhile, is crucial for the accuracy of yield estimation with remote-sensing techniques. Recent work recommends that combining different crop and non-crop parameters might significantly improve the accuracy of estimates (Choudhury et al., 2021; Jin et al., 2018). It is argued that vegetation indices alone may not capture all of the factors that affect crop yield, such as crop physiological and climate parameters of the crop (Cheng et al., 2022; Li et al., 2019). Some previous approaches have integrated various parameters for crop yield estimation models, such as soil moisture (Chakrabarti et al., 2014), evapotranspiration (Huang et al., 2015), solar radiation (Yu et al., 2020), or crop height (Li et al., 2011). Other studies combined hyperspectral, thermal and LIDAR UAV sensors and crop height to estimate corn, wheat and grass yields (Furukawa et al., 2020; García-Martínez et al., 2020). Few studies combined vegetation indices with various crop physiological and climate parameters of crops (Yu et al., 2020).

In addition to primary vegetation indices and plant height indicators, another higher-order indicator is solar radiation, which has a direct impact on crop growth, photosynthesis process, temperature regulation and water management (Cheng et al., 2022; Zhang et al., 2022). Only few studies (Holzman et al., 2018; Ihuoma & Madramootoo, 2020) so far have been combining vegetation indices and solar radiation indicators to improve accuracy of yield estimation model and reduce modelling errors. Their yield estimation models have demonstrated that the solar radiation indicator is extremely important to account for the energy input available for the crop to carry out photosynthesis efficiently. Existing models were, however, limited to solar radiation based on empirical formulas or satellite data, which needs to be calibrated for each study area. Moreover, these studies utilised MODIS (250 × 250 m) and GLDAS (>27.5 km) satellite data, which are not suitable for farm and field scales due to their relatively low resolution. Furthermore, the quality of pictures depends on suitable weather conditions at the requested point of time.

The majority of these studies, however, is based on expensive and complex sensors such as hyperspectral sensors, thermal sensors and Light Detection And Ranging (LIDAR) sensors (Furukawa et al., 2020; Tao et al., 2020; Westoby et al., 2012). Also, the few studies that integrated vegetation indices with crop height and other parameters (Choudhury et al., 2021; Panday et al., 2020; Tao et al., 2020) were all relying on named high-cost sensors or weather station data. UAVs equipped with simple multispectral cameras, also allow to collect spectral imagery, crop height, solar radiation and other crop physiological parameters of the crop (Choudhury et al., 2021), meanwhile at significantly lower cost. Low-cost crop yield estimation is particularly pertinent to low-income farmers, who are typically located in areas heavily affected by climate change anyway (Barnwal et al., 2013), but at the same time have low capital endowment and/or technical capacities for high-technology crop prediction.

This study’s aim is to test the options to increase wheat yield estimation efficiency exclusively with low-cost multispectral UAV imagery along two specific objectives: a) combining vegetation index with solar radiation and crop height indicators from multispectral UAV imagery to evaluate the accuracy of wheat estimations and b) comparison of wheat grain filling and early maturity stages to identify the period that provides a more accurate yield estimation with multispectral UAV imagery. With this approach, we contribute to the growing body of literature arguing for a multipronged and applicable approaches towards yield estimation instead of just focusing on expensive and complex equipment (Possoch et al., 2016). In practical terms, our findings will benefit farmers, agro-companies, insurance companies, as well as other agricultural stakeholders in low- and middle-income countries who cannot rely on high-cost technology for crop estimation. An UAV-based estimation approach might allow timely managerial decision to improve wheat yields also in low-income countries, where smallholder gain access to this technology via affordable drone service providers (Cook, 2017; Green et al., 2019). Ultimately, such improvements may contribute not only to attain stability of rural incomes as well as combating rising food prices and food scarcity.

Materials and methods

Study area

On-farm data collection was carried out in the southern part of Germany (Figure 1). The geographical location of the research areas is at the area of latitude 49°58’ North and longitude 11°10’ East at an elevation of about 570 m above the mean sea level.¹ This region has productive soil with a loamy texture, which prompts intensive use for agriculture (Flessa et al., 2002). Between 1991 and 2020, the annual mean temperature at the location was at 9.4°C (national average 9.3°), and the annual cumulative precipitation at about 740 mm (national average 808 mm) (Wetterdienst, 2022).

Figure 1. Three-dimensional view of one of the study fields.

As shown in Table 1, data were collected for three study fields in total, each between 3 and 6 hectares, which were very similar in terms of soil type and plot quality. Plots 1 and 2 were used to evaluate the accuracy of vegetation indices, solar radiation and crop height characteristics for wheat yield estimations in 2020. Plot 3 was used to determine the out-of-sample validity of our estimations.

Table 1. Plot characteristics.

Plot	Total acreage	Crop	Soil type	Plot quality index^a
Plot 1	3.20 ha	Wheat	Silty clay	34–45
Plot 2	6.15 ha	Wheat	Silty clay	44–48
Plot 3	5.2 ha	Wheat	Silty clay	36–44

^aPlot quality according to German soil quality indexing system, which is based on factors like climatic conditions, slope, shadow from forest etc. and ranges from 0–100, 100 being the highest quality soil in Germany.

General workflow

In this research, we examined the effectiveness of combination crop height, vegetation indices and solar radiation to estimate wheat yield based on multiple linear regression and quantile regression for 20 May (grain filling) and 16 June (early maturity) in 2020. We validated these results with data from a third plot in close vicinity. Analysis and processing of the collected data for the purposes of creating an ortho-mosaics, digital elevation models (Lobell et al., 2010) and digital surface models (DSM) were carried out with the software Agisoft PhotoScan Pro. The ready-to-use DEM and DSM data are calculated using Equation 1(1) $DEM-DSM=CSM$(1) to determine the wheat height parameter. The solar radiation parameter for 20 May (grain filling) and 16 June (early maturity) was calculated within the software ArcGIS 10.8 based on the DSM. The six VIs were selected based on their potential for wheat yield estimation determined in the previous research. In general, we examined the suitability of the Normalized Difference Vegetation Index (NDVI), Normalized Difference Red Edge index (NDRE), Chlorophyll red-edge Index (CIred-edge), Simplified Canopy Chlorophyll Content Index (SCCCI), Green Chlorophyll Vegetation Index GCVI and Enhanced Vegetation Index (EVI2) in combination with crop height (CH) and solar radiation parameters. To examine the influence of integrating VIs, CHs and solar radiation for crop yield estimation, we applied a multiple linear regression (MLR) and quantile regression. In the final step, we predicted wheat yields based on the best predictor variables and compared them with actual wheat yields.

The aforementioned steps are illustrated by Figure 2 and described more precisely in the subsections.

Figure 2. Methodological steps for wheat yield estimation based on UAV.

Remote sensing data collection

A multi-rotor unmanned aerial vehicle (UAV) was used to acquire ultrahigh-resolution multispectral images at an altitude of 50 m above the crop canopy with an image resolution of 4.5 cm, and a flight speed of 4 m/s (Figure 3). The maximum UAV flight operation time is 55 min. To configure and control the UAV flight trajectory, we used the software Mission Planner.

Figure 3. Multi rotor UAV.

The UAV is equipped with a five-band multi-spectral camera, namely RedEdge-MX by the company MicaSense (Figure 4). The RedEdge-MX camera outperforms conventional RBG cameras by including five multispectral bands (Table 2). The RedEdge-MX camera is set on a gimbal in order to mitigate the wind effect and thereby capture high-quality 3D images in TIFF format.

Figure 4. RedEdge-MX camera.

Table 2. Spectral wavelength information of the RedEdge-MX camera.

Band	Band Name	Wavelength (nm)	Bandwidth (nm)
1	Blue	475 nm	32 nm
2	Green	560 nm	27 nm
3	Red	668 nm	14 nm
4	Red-Edge	717 nm	12 nm
5	NIR	842 nm	57 nm

Wheat phenology consists of eight main growth phases: (1) germination and emergence; (2) early vegetative; (3) tilling; (4) flag leaf initiation; (5) booting; (6) anthesis; (7) grain filling and (8) maturity (Kaur et al., 2015; Siebert & Ewert, 2012). The optimum time period for using vegetation indices to estimate wheat, maize and soybean yields is when the crops are at the peak of their greenness level, in anthesis and grain filling phases (Guan et al., 2019; Kanning et al., 2018; Marino & Alvino, 2020). Furthermore, some studies indicated that the optimum time period for yield estimation based on crop height is the peak of the greenness level (anthesis and grain filling) as well as the maturity phase (Kanning et al., 2018; Song et al., 2020). Based on these literature findings, the local crop calendar, and farmers’ assessment of local growth stages, we selected 20 May (grain filling) and 16 June (early maturity) for UAV data collection in 2020. The collected imagery data comprises between 435 and 1000 filed pictures per plot.

Ground-level yield data

In order to validate our modelled data, we required reliable ground yield data. As shown by previous studies, self-reported crop yield estimates can be of low accuracy, leading to potentially biased results for training remote sensing algorithms (Paliwal & Jain, 2020). Instead, we relied on ground-level sensor data collected with a combine harvester, which is equipped with a grain yield monitor and real-time kinematic (RTK) global navigation satellite systems (GNSS) to record the dry crop weight and measure the harvested mass flow, moisture content, and geographic location (Song et al., 2020). The ground yield data from the harvester (example in Figure 5) were processed by the software ArcGIS (version 10.8). The harvester collected wheat yields from three wheat fields (Table 1), one of which was only used for out-of-sample validation.

Figure 5. Wheat yield map compiled from combine harvester machine.

In the final step, we generated a grid with pixel size 6.5 × 6.5 m a grid using ArcGIS (create fishnet tool); the total number of pixel observation was n = 1722.

Extraction of crop surface model for calculation plant height

UAV-obtained crop height data have the advantages of cost-effectiveness and flexibility, as compared to data obtained by hand measurement or a harvester during crop-cutting. The only scientific concern might be the accuracy of obtained crop height parameters: Two previous studies (Yu et al., 2020) measured the accuracy of a crop surface model (CSM) obtained with UAV imagery as compared to actual ground crop height, finding a very low estimation error in terms of mean crop height (<10%) as well as very high correlation between the two indicators (R² > 0.95). To determine the wheat height parameter, we obtained a digital elevation model and a digital surface model (DSM) with our UAV. For this research, the digital elevation model represents the elevation of soil at the pre-plant stage, e.g. soil under the crops. DEM data were collected twice, once at the pre-planting stage and once at the post-harvest stage to control for biases created by soil preparation. Meanwhile, we could observe no changes in elevation between both observation times. In the next step, we obtained a digital surface model (DSM), which is defined as the difference between the DEM and crop elevation, which was also collected with our UAV in May and June (Song et al., 2020). The DEM and DSM from the surveyed plots were extracted using Agisoft Metashape software. Agisoft Metashape is commonly used to generate DEM and DSM models, including mesh models, from a set of overlapping photographs (Khalaf & Hameed, 2020; Lastilla et al., 2021). When constructing a mesh using a dense point cloud, interpolation plays a part to create a continuous and smooth mesh surface. By utilizing triangular interpolation methods, the software can create high-quality meshes even in areas with missing data or occlusions, ensuring a more comprehensive and realistic representation of the captured area (Khalaf & Hameed, 2020; Lastilla et al., 2021). Finally, to extract real wheat plant height, we calculated a crop surface model (CSM) following Formula 1:

(1) $DEM-DSM=CSM$(1)

This data was then used to compute crop height from DEM and DSM. Within our sample plots, crop height was found to be between 0.35–1 m on 20 May, and 0.5–1.15 m on 16 June, averaging 0.93 in May and 0.87 in June, as also illustrated by the box plot in Figure 6.

Figure 6. Crop height for the wheat growing season.

Solar radiation

Solar radiation, defined as the amount of radiant energy for a given location, is another indicator used for crop estimation due to its importance for photosynthesis, hence plant growth and biomass accumulation and therefore directly affecting crop yields (Zhang et al., 2022). However, so far it has mostly been employed using solar radiation data from weather stations (Fuentes et al., 2020; Zhang et al., 2022). The solar radiation parameter from the surveyed plots was obtained from a digital surface model and calculated for 20 May and 16 June 2020 using the Solar Radiation tool of ArcGIS’s Spatial Analyst Toolbox. This tool calculates the solar radiation with elevation, slope, orientation and atmospheric transmission as most relevant inputs. The initial digital surface model was obtained from our 3D UAV model for the study area.

Vegetation indices

Finally, six different VIs were entering the model. Namely, we tested NDVI, NDRE, CIred-edge, SCCCI, GCVI and EVI2 (Table 3). All these parameters can be calculated with a comparatively affordable 5-band camera like the RedEdge-MX. The concentration of chlorophyll in the crops directly reflects crop health, showing a high correlation with the yield in the early ripening stage (Mamrutha et al., 2017; Tucker, 1979).

Table 3. Selection vegetation indices equation and their application using UAV bands.

Vegetation Index	Equation	Reference
Normalized difference vegetation index (NDVI)	$NDVI=\frac{P_{NIR}-P_{RED}}{P_{NIR}+P_{RED}}$	C. J. J. R. s. o. E. Tucker (1979)
Normalized difference red edge index (NDRE)	$NDRE=\frac{P_{NIR}-P_{RED-Edge}}{P_{NIR}+P_{RED-Edge}}$	Gitelson and Merzlyak (1994)
Chlorophyll red-edge index (CIred-edge)	$C{I}_{red-edge}=\frac{P_{NIR}}{P_{RED-Edge}}-1$	Gitelson et al. (2003)
Simplified canopy chlorophyllcontent index (SCCCI)	$SCCCI=\frac{\frac{P_{NIR}-P_{RED-Edge}}{P_{NIR}+P_{RED-Edge}}}{\frac{P_{NIR}-P_{RED}}{P_{NIR}+P_{RED}}}$	Barnes et al. (2000)
Green chlorophyll vegetation index (GCVI)	$GCVI=\frac{P_{NIR}-P_{Green}}{P_{NIR}+P_{Green}}$	Gitelson et al. (2003)
Enhanced vegetation index 2 (EVI2)	$EVI2=2.4\frac{P_{NIR}-P_{RED}}{P_{NIR}+P_{RED}+1}$	Jiang et al. (2008)

Regression model

We applied several linear regression models to examine the influence of integrating VIs, solar radiation and CHs for crop yield estimation as implemented in previous studies (Bendig et al., 2015; Thenkabail et al., 2000). For multi-indicator estimations, our multi-linear-regression (MLS) model can be written as:

(2) $y_{\hspace{0.17em}}={\beta }_0+{\beta }_1{x}_{CH}+{\beta }_2{x}_{VI}+{\beta }_3{x}_{SR}+\epsilon$(2)

where y is the ground actual wheat yield, and x_VI, x_SR, and x_CH are the explanatory variables with their model coefficients β. In our case, x_VI, x_SR, and x_CH are the vegetation indices (VIs), solar radiation (SR) and crop height (CH), respectively, and $\epsilon$ is defined as the error term.

Furthermore, we also examined the suitability of the VIs, solar radiation and CHs variables separately based on a simple linear regression model, written as:

(3) $y_{\hspace{0.17em}}={\beta }_0+{\beta }_{1}x+\epsilon$(3)

where y is the ground actual wheat yield, x represents the selected VI, CH or radiation indicator and $\epsilon$ is defined as the error term.

All regression analyses were conducted with Stata 16.0 software. The model performance was evaluated along the Root Mean Square Error (RMSE), the R squared (R²) and the adjusted R² (Table 4).

Table 4. Definitions and equations of precision evaluation metrics used in this study.

Name	Equation
R squared (R^2)	$R^2=1-\sum i=1n{y}_i-{\widehat{y}}_{i}2{\displaystyle {\sum }_{i=1}^n{\left(y_i-\overline{y_i}\right)}^2{ }^{\prime }}$
Root Mean Squared Error (RMSE)	$RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(y_i-{\hat{y}}_i\right)}^2}$

Quantile regression was utilized to analyze different ranges of crop yield at different minimum, median and maximum quantiles. Usually, studies on crop yield estimation mainly focused on the mean and variance parameters of crop yields. The disadvantage of this approach, however, is that estimations are extremely sensitive to outliers and the requirement of independence of the residuals is not satisfied (Mishra & Moss, 2013; Radhakrishna & Toutenburg, 1995). Also, crop yield estimation based on ordinary least squares (traditional regression methods) may not give adequate results due to extreme values or non-normal distribution in climate data (Murungweni et al., 2020; Tareghian & Rasmussen, 2013). Furthermore, the distribution of errors in the lower or upper tails of crop yield distribution is significant for precise crop estimation, especially when applied for drought early warning systems or agricultural insurance.

Quantile regression helps to overcome the possible asymmetry in impacts across each quantile of conditional crop yields over certain pivotal variables (Barnwal et al., 2013). Another huge benefit of quantile regression is its robustness to outliers (Koenker & Bassett, 1978; Mishra & Moss, 2013), which is an important factor in estimating crop yields. Moreover, quantile regression can enable the detection of all nuanced views of the stochastic relationship among parameters, and therefore a more detailed empirical analysis.

Our quantile regression equation is given below:

(4) $Q_{\tau }\left(y_{\hspace{0.17em}}\right)={\beta }_0\left(\tau \right)+{\beta }_1\left(\tau \right)x_{CH}+{\beta }_2\left(\tau \right)x_{VI}+{\beta }_3\left(\tau \right)x_{SR}+\epsilon \left(\tau \right)$(4)

Here, $y_{\hspace{0.17em}}$ is the conditional value of the actual wheat yield, whereas β₁ and β₂ are the estimation parameters, and τ indicates the quantile (i.e. τ = 0.25, 0.5, 0.75 for the minimum, median, maximum yield). x_CH, x_SR and x_VI are the value of the independent variables VIs, SR and CH, and $\epsilon$ are the error terms (Murungweni et al., 2020; Vilar et al., 2015).

We apply quantile regression to determine the prediction power of vegetation indices, solar radiation and crop height for wheat yields at different quantiles (0.25, 0.50, 0.75), e.g. lower, median and higher yield. In line with Koenker & Machado (1999), the goodness-of-fit of the regression model for each quantile (0.25, 0.50 and 0.75) is demonstrated by the value of Pseudo R². Moreover, we generated model residual plots to evaluate how effectively the model predicts crop yield (Elavarasan et al., 2020).

Results

Wheat yield estimation based on linear regression models

Table 5 displays a number of simple linear regression models measuring the correlation between a range of separate vegetation indices and ground-level yield data. We observe a moderate but significant correlation between actual yields and VI on 20 May, in particular for the NDRE index (R² = 0.60, adj R² = 0.60, RMSE = 0.68). On 16 June, we only measured a moderate correlation between GCVI and actual yield (Figure 7, R² = 0.50 and RMSE = 0.77). Individual crop height or solar radiation variables meanwhile had lower accuracy (see Supplementary Table S1, R² < 0.36, adj R² < 0.35, RMSE > 0.88) than VIs for estimating wheat yields. Furthermore, our analysis reveals a higher accuracy of indices using the near-infrared and the red edge bands (NDRE, CIred-edge, SCCCI) as compared to indices using a combination of the near-infrared band, the red band and the green band (NDVI, EVI2, GCVI) (see Supplementary Table S1). The results of the simple linear regression therefore suggest a high sensitivity of a red edge band (NDRE) to wheat health changes during the maturation period. Furthermore, NDRE revealed a strong correlation at the peak of greenness level (before early maturity) compared to the remaining vegetation indices (Figure 7).

Figure 7. Results of the multiple and linear regression models a) R² and b) RMSE.

Table 5. Validation results MLR between actual wheat yield data and combined NDRE_May20, SR_June16 with CH_June16.

yield_1	Coef.	St.Err.	t-value	p-value	[95% Conf	Interval]	Sig
ndre_20	17.914	0.538	33.32	0	16.858	18.969	***
csm_16	1.242	0.068	18.39	0	1.109	1.374	***
sr_16	−1.994	0.618	−3.22	0.001	−3.208	−0.78	***
Constant	6.389	3.269	1.95	0.051	−0.028	12.806	*
Mean dependent var	5.736		SD dependent var			1.069
R-squared	0.75		Number of obs			1722.000
Adj R-squared	0.74		Root MSE			0.53
F-test	715.492		Prob > F			0.000
Akaike crit. Aich et al. (2017)	1152.935		Bayesian crit. (BIC)			1171.264

*** p < 0.01, ** p < 0.05, * p < 0.1.

In the second step, we measured the relationship between actual wheat yield data and the combination of VIs, solar radiation and CHs along a multiple linear regression model displayed in equation (2)(2) $y_{\hspace{0.17em}}={\beta }_0+{\beta }_1{x}_{CH}+{\beta }_2{x}_{VI}+{\beta }_3{x}_{SR}+\epsilon$(2) . The results of which are displayed in Figure 7 for 20 May and 16 June. All combinations were statistically significantly correlated. The highest correlation level was found for a combination of CH, solar radiation and the NDRE index (R² = 0.75 and RMSE = 0.53) for CH and solar radiation on 16 June and the NDRE index on 20 May. A similar precision was found for a combination of CH, solar radiation and the CIred-edge index (R² = 0.72 and RMSE = 0.56) where CH and solar radiation were measured on 16 June and the NDRE index was measured on 20 May. All other combinations of VI, solar radiation and CH parameters for both time observations revealed a medium correlation (Figure 7, GCVI and SCCCI) to low correlation (Figure 7, NDVI and EVI2). The comparative analysis of the VI regressions with actual wheat yield undertaken on different days reveals a higher correlation on 20 May, when the wheat canopy reflectance has the highest greenness (Figure 7). A lower explanatory power of VIs for wheat yield was observed on 16 June, when wheat is at an early maturity stage, in other words starting to turn yellow (Figure 7). In return, wheat CH and solar radiation parameters were more strongly correlated with wheat yields on 16 June than in 20 May, a time when length growth peaks (Figure 7).

Combining vegetation indices with crop height and solar radiation significantly improves the accuracy of the estimation model (Figure 7) as compared to utilizing them individually (Figure 7 and Supplementary Table S1), which is in line other findings (Cheng et al., 2022; Choudhury et al., 2021). Among the combinations between vegetation indices with other parameters, the highest performance was found for NDRE (Figure 7, R² = 0.75 and RMSE = 0.53) and CIred-edge (Figure 7, R² = 0.72 and RMSE = 0.56).

In summary, the highest performance was found for the estimation combining the NDRE index from 20 May with crop height and solar radiation characteristics from 16 June. In Table 5, we provide comprehensive regression results for this specific model. For the other models, detailed regression outputs have been shifted to the Supplementary Table S2.

Wheat yield estimation based on quantile regression models

Generally, we find a statistically significant relationship between wheat yield and crop height/solar radiation/vegetation indices for all quantiles. Meanwhile, the Pseudo R-squared (pR²) in all quantiles presented a non-high pR² coefficient, except for the NDRE index combination with crop height characteristics and solar radiation for 20 May and 16 June (see Supplementary Table S3). We also observe the strongest relationship between actual wheat yield and the combination of NDRE index with crop height and solar radiation for the lower yield quantile (Figure 8 and Table 6). The median and high yield quantile demonstrated a lower relationship with VIs than the low yield quantile (Figure 8). Overall, the highest correlation level was found for the combination of NDRE index from 20 May and CH and solar radiation from 16 June for wheat yields in the lower quantile (R² = 0.57). In summary, quantile regression results confirmed that the combined VIs, solar radiation and wheat CH parameters can predict wheat yields more accurately than predictions with one indicator separately.

Figure 8. Results of the quantile regression model.

Table 6. Quantile regression result between actual wheat yield data and combined NDRE_May20, SR_June16 with CH_June16.

yield	Coef.	St.Err.	t-value	p-value	[95% Conf	Interval]	Sig
ndre_20	18.754	0.782	23.99	0	17.22	20.289	***
csm_16	1.131	0.086	13.21	0	0.963	1.3	***
sr_16	−1.009	0.602	−1.68	0.094	−2.19	0.173	*
Constant	−2.023	3.079	−0.66	0.511	−8.068	4.022
ndre_20	17.515	0.4	43.82	0	16.73	18.3	***
csm_16	1.075	0.099	10.89	0	0.881	1.269	***
sr_16	−2.822	0.784	−3.60	0	−4.361	−1.283	***
Constant	8.93	4.346	2.05	0.04	0.398	17.463	**
ndre_20	16.393	0.512	32.02	0	15.388	17.398	***
csm_16	1.056	0.065	16.18	0	0.928	1.184	***
sr_16	−2.638	0.636	−4.14	0	−3.887	−1.388	***
Constant	8.804	3.461	2.54	0.011	2.009	15.599	**
Mean dependent var	5.736		SD dependent var			1.069
Pseudo r-squared 0.25	0.578		Number of obs			1722.000
Pseudo r-squared 0.50	0.544		Number of obs			1722.000
Pseudo r-squared 0.75	0.556		Number of obs			1722.000

*** p < 0.01, ** p < 0.05, * p < 0.1.

To provide a more precise picture, we additionally provide quantile plots. The magnitude of the quantile line demonstrates non-linear oscillations at lower quantiles from 0.1 to 0.2 (Figure 9). For other indicator combinations, we refer to Supplementary Figure S1. All quantile regression plots demonstrated that quantile regression is more appropriate compared to multiple linear models, especially in the lower quantiles of the yield (see Supplementary Figure S1). Meanwhile, linear regression models provide adequate estimates mainly beyond the second decile.

Figure 9. Estimated coefficients with 95% confidence intervals for NDRE_May20 , SR_June16 and CH_June16.

Out-of-sample validation

We also conducted an estimation of wheat yields for nearly plot for 2021 to test the validity of our estimation outside this specific sample and compared the estimated with true ground yield data. The results of this comparison are provided in Figure 10. Figure 10(a) demonstrates that the predicted yields for each pixel closely reflects the actual wheat yields. Figure 10(b) is a scatter plot of actual wheat yields and predicted yields based on the combination of NDRE index, solar radiation and CH parameters.

Figure 10. a) comparison of model-predicted wheat yield with actual yield; b) scatter plot of predicted wheat yield with actual yield.

Figure 10(a) demonstrates that the wheat yield prediction line does not present a serious inflated or undervalued estimation situation and the wheat prediction line has approximately 10–15% deviation from the actual yield line. In Figure 10(b), the R² coefficient is indicated on the left, which represents the proportion of the deviation of the predicted wheat yield from the actual yield. Here, the R-square ranges at 0.80, which means that the predicted and actual wheat yield values match by at least 0.80% (Figure 10(b)).

Figure 11 shows a scatter plot of the model residuals based on a combination of NDRE index and CH parameters with actual wheat yields. Ignoring a few outliers at the lower left corner of the plot, we can see that all model residuals are randomly distributed along a horizontal line, demonstrating adequate model performance (Sullivan, 2018). We might also have seen increasing or decreasing variation in the residuals’ heteroscedasticity.

Figure 11. MLR model residuals plot.

Discussion

In our results, we firstly presented the isolated performances of separate indices in wheat yield predictions. Among the tested indices, the NDRE performed best in predicting wheat yields. The high precision of the NDRE index can be related to its reliance on the reflectance of red-edge wave lengths, which proves to be more responsive to variations in leaf chlorophyll content than the red reflectance utilized by NDVI and EVI2. This heightened sensitivity of the NDRE index makes it better at detecting changes in crop health, a property that is directly linked to wheat yields, as highlighted in studies by Boiarskii and Hasegawa (2019) or Zeng et al. (2022). On the other hand, traditional indices like NDVI and EVI2 may respond less to rapid increases in above-ground biomass of wheat at the grain-filling stage of wheat growth. Several recent studies (Guan et al., 2019; Kanning et al., 2018; Song et al., 2020) demonstrated that the SCCCI and CIred-edge indices provide more accurate estimates of yield than other vegetation indices. In our study, however, these indices do not outperform the NDRE index in terms of accuracy, which may be related to the fact that these previous studies utilized hyperspectral and LIDAR sensors. This observation is an important hint that the vegetation index performance is not only dependent on location and data resolution but also on the type of sensors used for data collection.

Integrating NDRE index values, solar radiation and the wheat height parameter results in an improvement of up to 15–20% in predicting wheat accuracy compared to utilizing either of the VIs individually. Each of these indicators provides unique insights into various information related to crop health, energy availability, and physical development. While NDRE measures plant health, solar radiation considers energy accessibility and wheat height reflects the physical development of the plant. All of these measures provide a holistic view of crop health and growth when used conjointly, which explains the higher performance of the combined estimation.

A comparison of wheat growing stages revealed that solar radiation and crop height characteristics provided a more accurate estimate of crop yield at the early maturity stage (16 June), while vegetation indices provided a more accurate estimate at the grain filling stage (20 May). The wheat crop is nearing harvest in the early maturity stage (Hassan et al., 2019; Holman et al., 2016; Yuan et al., 2018), hence, its growth does not much change until being harvested. Furthermore, during this stage efficient photosynthesis can be observed. The high accuracy of vegetation indices at the grain filling stage is related to the fact that wheat has the highest green leaf reflectance at this stage. Moreover, at these stages, any stressors or problems affecting crop health can be easily detected.

One of our main motivations of this study was to test if shortcomings of low-cost UAV sensors can be evened out by a multi-indicator estimation. In facts, we obtained equal or higher precision of estimates (R² = 0.75, Adj. R² = 0.74 and RMSE = 0.53) than previous studies (R² > 0.70, Adj R² = 0.68 and RMSE = 0.60) that used high-cost hyperspectral or sensor data (Choudhury et al., 2021; Tao et al., 2020). Hence, we can confirm that combining the NDRE index value, solar radiation and CH parameters is an effective remote-sensing method to increase the accuracy of wheat yield estimates on a field scale even under low technology input. Thus, our approach allows the implementation of UAV-based crop estimation also for smaller family farms in low-income countries via service providers. It should also be pointed out that we achieved these results despite not working on controlled-environment university or laboratory fields but on real farm plots (Belton et al., 2019; Tao et al., 2020).

A number of limitations of the proposed approach must be recognized. Firstly, the research was conducted for one year and two vegetation times. While we tested the robustness within a second year, we would require multiple years and plots for further robustness tests. Secondly, the inclusion of additional indicators affecting yields, such as soil quality, temperature, precipitation, or fertilizer would allow incorporating differences in management practices, which should be considered when conducting analysis in large areas with heterogeneous management, however would once again render the estimation more complex. Another main prospect for future research is to combine long-term data with panel data regression, as well as incorporating a crop-cutting approaches to further confirm the validity of obtain crop height parameters.

Conclusion

In this study, we tested the suitability of integrated various physical and spectral information of wheat based on low cost multispectral images obtained from the UAV data images for wheat yield estimation. Combining the NDRE index with a wheat height parameter and solar radiation can significantly improve model accuracy, minimize the weaknesses of individual indicators, and decrease RMSE (Figure 7 NDRE, SR and CH RMSE < 0.53 and NDRE RMSE > 0.68) by 10–15%. Our three main findings are as follows: Firstly, multiple-indicator yield estimates in our case were more accurate than a single indicator, while the combination of the NDRE index with the wheat height parameter and solar radiation performed the best. Secondly, the growth time of wheat at grain filling is more appropriate for estimating yields based on vegetative indices, but the growth time of early maturing wheat is better using crop height and solar radiation indicators. Finally, utilizing low-cost sensors to estimate wheat yield demonstrated strong potential in terms of accuracy of yield predictions, also compared to high-cost benchmark estimation.

Our approach could contribute to improved decision-making in agricultural production and management, as well as climate change mitigation and adaptation. The UAV-based approach can empower farmers with valuable information on monitoring crop growth for optimizing yield potential and enhancing agricultural as well as land management practices. Moreover, estimations would help farmers to identify specific sites that require attention and implement targeted management practices as well as minimize potential damage to the yield.

Acknowledgments

The authors express their gratitude to the team members of the project “Increasing Climate Resilience via Agricultural Insurance – Innovation Transfer for Sustainable Rural Development in Central Asia (KlimALEZ)”, implemented by the Leibniz-Institute for Agricultural Development in Transformation Economies (IAMO).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/22797254.2023.2294121.

Additional information

Funding

This research was supported by the German Federal Ministry of Education and Research (BMBF) [01LZ1705A] and WGL Forschungsverbünde [D0093/3].

Notes

¹ For data security reasons (private farms) we only gave the approximate location of the study area.