Introduction

Spatial variability in crop yields on arable fields

Arable fields are characterized by a more or less pronounced spatial variability in crop yields (Schuster et al., 2023). One of the main causes of spatial variability in plant growth and biomass production is heterogeneous soil properties (Godwin et al., 2003; Hatfield, 2000; Mittermayer et al., 2021, 2022), such as soil texture, soil organic carbon content (SOC), total nitrogen content (TN), macro- and micronutrient content, and available water capacity (López-Lozano et al., 2010; Servadio et al., 2017). The spatial variation in soil properties is influenced by pedogenesis, topography and soil erosion (Gregory et al., 2015; Raimondi et al., 2010). Cultivation practices can also contribute substantially to yield variability (Ngoune & Shelton, 2020), e.g., through soil compaction by heavy machinery (Horn & Fleige, 2003; Shaheb et al., 2021) or management practices (fertilizer and pesticide application, weed control). Natural field boundaries through adjacent forest or tree strips and agroforestry systems also influence yield variability within fields (Karlson et al., 2023; Pardon et al., 2018). In addition, yield variability is influenced by weather conditions. In years with drought stress, yield zones are more pronounced than in years with better rainfall distribution (Heil et al., 2023; Maestrini & Basso, 2018; Martinez-Feria & Basso, 2020). Some factors influencing yield variability are stable over a long-term period (e.g., soil texture), whereas others change from year to year.

Consideration of spatial yield variation is an important factor for precision agriculture and the key to optimizing crop production through more efficient use of inputs (fertilizers and pesticides) (Gebbers & Adamchuk, 2010; Mulla, 2013).

Although some farmers are currently considering the use of precision farming applications with variable, yield-dependent application rates, the use of these digital technologies rarely exceeds more than 20% of farms, as farmers are not yet convinced of the benefits of these methods (Lowenberg‐DeBoer & Erickson, 2019; Gabriel & Gandorfer, 2022). Considering that uniform nitrogen fertilization of arable fields is still a common practice, spatial variability in yields may lead to nitrogen oversupply in low-yield zones, which can cause high nitrogen losses, e.g., nitrate leaching (Mittermayer et al., 2021; Schuster et al., 2022). In general, the yield potential should be accounted for to optimally adapt crop management measures, such as fertilization and plant protection accordingly, which requires yield analysis methods that are as precise, cost-effective and widely available as possible.

Research needs

To delineate yield zones for precision farming applications, knowledge of the spatial variation of crop yield is crucial. In agricultural practice, digital combine harvester yield sensing systems are most commonly used for the determination of yield variability within fields, although raw data from combine harvesters have a large error potential (Bachmaier, 2010; Fulton et al., 2018; Kharel et al., 2019). Ground truth data on spatial variation in grain yield may be collected using nondigital methods, such as georeferenced plant sampling or plot harvesting with a plot combine harvester (Mittermayer et al., 2021; Spicker, 2017; Stettmer et al., 2022b). These methods are expensive and labor intensive and can therefore only be used for scientific studies.

Several studies have been conducted on spatial yield estimation based on multispectral measurements by remote and proximal sensing (Aranguren et al., 2020; Barmeier et al., 2017; Maidl et al., 2019). Based on reflectance data, vegetation indices (VI), such as the normalized difference vegetation index (NDVI), were found to be closely related to grain yield, plant biomass and nitrogen uptake at specific growth stages (Cabrera-Bosquet et al., 2011; Prabhakara et al., 2015). Various methods for spatial yield estimation have been developed using artificial intelligence (AI) (Ruan et al., 2022; van Klompenburg et al., 2020).

However, most methodical approaches for digital spatial yield analysis are not fully transparent and comprehensible because the algorithms used are not described in sufficient detail. In some cases, commercial suppliers use very complex soil process and plant growth models that are difficult for users to understand. Moreover, important agronomic parameters (e.g., crop type, growth stages and weather conditions) are neglected in some satellite-based yield estimation systems (Li et al., 2007). Furthermore, in some cases, the methods have not been sufficiently validated; thus, the accuracy of the delineated zones is questionable.

Independent validations showed that in some cases substantial deviations in the satellite-derived yields from measured yields can occure (yield differences by several Mg ha−1) (Mittermayer et al., 2021; Stettmer et al., 2022b). Therefore, management decisions should not be based on such large deviations between estimated yields and actual yields.

Most studies analyzing spatial yield variability focus on plant-based variables (e.g., grain yield, biomass, nitrogen uptake), while soil-related causes are neglected. Due to high costs, high-resolution soil mapping is not often conducted, and thus, important influencing factors that are partly responsible for the variability in yield and biomass formation of crops are not accounted for (Feng et al., 2022; Juhos et al., 2015).

Although there are numerous scientific studies on satellite-based yield analysis, as well as commercial service providers selling yield maps to farmers, there is a need to further develop satellite-based yield analysis methods, improve their accuracy and validate them under differentiated management conditions.

Study aims

This study describes a new approach to generate relative biomass potential maps (rel. BMP). In this context, relative biomass potential means the determination of multiyear stable yield zones within fields at a high spatial resolution, without information on the absolute yield level in the yield zones. For many applications, knowledge of relative yield differences is sufficient; estimating absolute yields requires more input data and more complicated models.

Data (time series of several observation dates) from Sentinel-2 satellites of the European Space Agency (ESA), the vegetation index NDVI and agronomic knowledge were used to determine the relative biomass potential.

Over a twenty-year research period at the Technical University of Munich, multispectral proximal sensor systems have been systematically used to determine the vegetation indices, and at what vegetation stages, that correlate best with yield, nitrogen uptake, and biomass formation (Mistele & Schmidhalter, 2010; Mistele et al., 2004; Schmidhalter et al., 2003a; Spicker, 2017; Sticksel et al., 2004). Analysis of the relationships between vegetation indices and agronomic parameters has shown the importance of considering the correct growth stages of plant stands as influencing factors (Erdle et al., 2011; Prey & Schmidhalter, 2019; Schmidhalter et al., 2003b)). Therefore, one aim of this study was to test whether the knowledge gained from previous proximal sensor-based studies, the relationships found between vegetation indices, and agronomic parameters can be transferred to satellite-based yield estimation methods. Satellite data from observation dates corresponding to specific growth stages were used to analyze crop-specific yield patterns. For winter wheat, for example, the growth stages jointing (GS 32), booting (GS 39) and flowering (GS 65) were considered (AHDB, 2023; Zadoks et al., 1974). To achieve a high accuracy of the biomass potential estimates, data from many observations and different crop types in the crop rotation were analyzed and combined. In this study, the crop species winter wheat (Triticum aestivum L.), winter barley (Hordeum vulgare L.), canola (Brassica napus L.), corn (Zea mays L.) and soybeans (Glycine max L. Merr.) were analyzed. Due to yearly changing weather conditions, site-specific yield patterns may vary from year to year, and yield maps can be year specific. Therefore, the rel. BMP maps were analyzed for years with different weather conditions (dry, wet, normal conditions), and whether a multiyear data analysis increases accuracy and validity was investigated.

An analysis of the relationship among yield patterns with spatial variability in soil property contents (e.g. SOC, TN, phosphorus (P), potassium (K), pH and soil texture) was part of this study. In addition, multispectral reflectance measurements were conducted with a tractor-mounted sensor system to obtain high-quality validation data using existing and validated yield algorithms, independent of the satellite data.

Based on the current state of knowledge, the following hypotheses were formulated:

  1. (1)

    Consideration of the crop type and crop-specific growth stages increases the accuracy in deriving yield zones compared to methods that do not include this information.

  2. (2)

    The yield patterns of the respective crops are similar if the correct crop-specific growth stages are selected.

  3. (3)

    When deriving yield zones, satellite data almost reach the same accuracy as data from tractor-mounted sensor systems.

  4. (4)

    The spatial distribution of crop-specific relative biomass yield potential is closely related to the spatial distribution of ground truth data (e.g., biomass yield, soil properties influencing yield).

Materials and methods

Site and weather conditions

The investigations were conducted on arable fields at three different sites. The methodology for compiling biomass potential maps was derived from fields at two research stations (Roggenstein (site A) and Dürnast (site B)) of the Technical University of Munich (Bavaria, Germany) and was validated on arable fields of farms in the Burghausen region (80 km east of Munich 48° 7′ 51″ N 12°44′ 5″ E, 450 m a.s.l., site C) under different management conditions. Research station Dürnast is located 30 km north of Munich (48° 24′ 3″ N 11°38′ 42″ E, 425 m a.s.l.) and research station Roggenstein is located 20 km west of Munich (48° 10′ 47″ N 11° 18′50″ E, 480 m a.s.l.).

Due to the high availability of yield data (determined with different systems) and many years of precision farming experience, the research stations were used to derive the relative BMP algorithm. Sites A and B are located in the tertiary hill country of southern Germany. The major soil types are cambisols (medium to deep brown soil of medium quality) (FAO, 2014). The 30-year mean annual precipitation was 888 mm (site A) and 813 mm (site B), and the temperature was 8.9 °C (Online Resource 1, Online Resource 2). Site C is located on the Alzplatte, which is characterized by a smooth, hilly landscape. The 30-year mean annual precipitation was 849 mm, and the mean annual temperature was 8.9 °C (Online Resource 3).

To derive and validate the satellite-based method for analyzing relative yield potential, multiyear data (2018–2022) were used to account for different weather conditions in the study regions. A dry and warm spring and a hot summer characterized 2018 and 2019. At all sites, the mean precipitation in 2018 and 2019 was 10 to 20% below average. In 2020, the spring was dry at all sites, while the summer was within the normal range with some heavy rainfall. The year 2021 was very wet with heavy rainfall. The temperatures were average. In 2022, the spring was dry, and the summer months were warmer with less precipitation (12 to 25%) than the 30-year average.

Farm management of the investigated fields

The investigated fields were cultivated in conventional arable farming systems. All fields were at least four hectares in size (above average size in the regions). The crop rotation of each field is shown in Table 1. In fields A1 (site A) and B1 (site B), only mineral fertilizer was applied. In fields C1 (site C) and C2 (site C), organic fertilizers have been applied for several years. The major crops of all farm sites were winter wheat, winter barley, corn, canola and soybeans (Online Resource 4).

Table 1 Study fields
Full size table

Remote sensing data

The algorithm was developed using Sentinel-2 MSI-Level 2A (MAJA Tiles) satellite data provided by the German Aerospace Center (DLR). The Sentinel-2A mission from the Copernicus program by the European Space Agency (ESA) provides satellite images with 13 spectral bands, four of them (red (665 nm), green (560 nm), blue (490 nm), and VNIR (842 nm)) at a 10 × 10 m resolution. The return rate of the two satellites of the Sentinel-2 mission is 5 days (Sentinel Online, 2023). The satellite data were preprocessed using the time-series-based MAJA processor, a multitemporal atmospheric correction and cloud screening algorithm developed by the DLR (German Aerospace Center, 2019; Hagolle et al., 2017).

Description of the relative biomass potential map algorithm

To develop the algorithm, satellite data from 2018 to 2022 were retrospectively acquired from DLR. To display within-field spatial variation, a 5-year time series of Sentinel-2A images was processed. As the vegetation index NDVI has already been demonstrated to estimate plant aboveground biomass (Kross et al., 2015; Perry et al., 2022), NDVI values were used as an indicator of the biomass growth potential on cropland. An essential aspect of the algorithm is the knowledge of the cultivated crop types and the dates of the characteristic development stages. This information was available in detail for the fields to derive the algorithm at the research stations, but also on the validation fields of farmers. An overview of the workflow for determining the rel. BMP map according to Hagn et al. (2023) is shown in Fig. 1. The work steps are described as follows:

Fig. 1
figure 1

Workflow of the rel. BMP (%) map algorithm, presented using the example of winter wheat

Full size image

Step (1): Selecting the satellite scenes according to the grown crops

Depending on the crops grown in the field, satellite scenes at characteristic growth stages are selected (Online Resource 5). The growth stages of booting (GS 30–32), jointing (GS 39) and flowering (GS 65) of winter wheat and winter barley correlate well with yield and biomass growth (Barmeier et al., 2017; Erdle et al., 2011; Maidl et al., 2019; Mistele & Schmidhalter, 2008a; Prey & Schmidhalter, 2019; Spicker, 2017). However, for corn and soybeans sown in split rows, growth stages have to be selected at which row closure has already been reached, for corn the growth stages GS 39 and GS 65 (Mistele & Schmidhalter, 2008b), and for soybeans the growth stages V5 (fifth trifoliate), R2 (full bloom) and R5 (beginning of seed filling) (Crusiol et al., 2016; Farias et al., 2023). For canola, the growth stages budding (GS 50), beginning of flowering (GS 60), where no yellow coloration occurs yet, and podding (GS 70), where the yellow coloration of the flowers has already faded again, are used (Spicker, 2017).

Step (2): Clipping of Sentinel-2A data to the field and checking for clouds

The satellite scenes are clipped to the field. Then, for each satellite scene date, the NIR and red bands are used to check whether the area to be analyzed was covered by clouds. If so, the image is rejected.

Step (3): Calculation of the NDVI

The vegetation index NDVI is calculated for each 10 × 10 m raster cell of each satellite scene according to the growth stages, and the mean NDVI of the whole field of each satellite scene is calculated. Steps (1) to (3) are repeated for all growth stages of the cultivated crop mentioned in Step (1).

$$NDVI = \frac{842 nm - 665 nm}{\ 842 nm + 665 nm}$$
(1)

Step (4): Calculation of the relative NDVI (rel. NDVI)

The rel. NDVI (%) is calculated by dividing the NDVI of each 10 × 10 m raster cell of each satellite scene by the mean NDVI (\(\overline{NDVI }\)) of the entire field of each satellite scene and multiplying it by 100. Meaning that in total three rel. NDVI maps are generated per year for each crop (exept for corn) of the crop rotation.

$$rel. NDVI \left( \% \right) = \frac{NDVI}{{\overline{ NDVI} }}*100$$
(2)

Step (5): Calculation of the singleyear rel. BMP (%)

The singleyear rel. BMP is determined by calculating the arithmetic mean of all relative NDVI values of the raster cells of the rel. NDVI maps generated in Step (4). In the case of wheat, the rel. NDVI values of the respective grid cells of the rel. NDVI maps of GS 32, GS 39 and GS 65 are summarized and divided by their number. Steps (4) and (5) are repeated for all crops of the crop rotation, or at least for a period of five years, resulting in five singleyear rel. BMP maps.

$${\text{Singleyear rel}}.{\text{ BMP }} = \mathop \sum \limits_{i = 1}^{n} rel. NDVI_{i} = { }\frac{{rel. NDVI_{GS32} + rel. NDVI_{GS39} + rel. NDVI_{GS65} }}{3}$$
(3)

Step (6): Calculation of the multiyear rel. BMP

The multiyear rel. BMP (%) is determined by calculating the the arithmetic mean of all singleyear rel. BMP maps. For a five year crop rotation, five singleyear rel. BMP are created in Step (5), meaning that the total of all singleyear rel. BMP are divided by five.

$$Multiyear rel. BMP = \mathop \sum \limits_{i = 1}^{n} rel. BMP_{i} = { }\frac{{rel. BMP_{1} + rel. BMP_{2} + rel. BMP_{3} + rel. BMP_{n } }}{n}$$
(4)

Step (7): Meaning and interpretation of the rel. BMP map

The map of relative biomass potential reflects the multiyear site-specific yield potential of a crop rotation. By accounting for different crops and observation dates, greater precision should be achieved than with a single analysis of a specific crop. The primary objective of the relative BMP maps is the identification of different yield zones within arable fields to identify high- and low-yield zones for site-specific management.

Methods of determining grain yield

At field A1 and field B1, the yield was measured by combine harvesters with an integrated yield sensing system (Claas Lexion 5500 and New Holland CX 790). The grain yield was determined by three main components (GPS-position with Real Time Kinematic, determination of the harvested area, determination of the grain moisture and the grain yield by means of a moisture sensor and volume flow measurement sensor (Noack, 2006).

At fields C1 and C2 the yield was derived based on the spectral reflectance measurements from the tractor-mounted sensor system in combination with a yield algorithm according to Maidl et al. (2019).

As a nondigital method for grain yield determination, 50 georeferenced biomass samples were taken before grain harvest at field A1. Eight 2 m long winter wheat rows were cut close to the ground with hand shears. The samples were threshed with a stationary combine thresher (Wintersteiger, 2023). After drying the grain at 60 °C, the dry matter (DM) content and the yield (t ha−1) at 86% DM were determined.

Methods of determining spatial reflectance data by proximal sensing

Spectral reflectance measurements during the flowering of winter wheat were carried out with a tractor-mounted multispectral sensor system (TEC5 2010) in 2018, 2021 and 2022. The system is equipped with two multispectral reflectance sensors (360 nm to 900 nm). Approximately every second reflectance measurements are conducted. Natural influences of solar radiation are considered in the data output by the implementation of a reflectance reference module. Based on the reflectance measurements, the vegetation index REIP (red edge inflection point) was calculated (Guyot et al., 1988). Since various studies have shown that REIP is closely related to the aboveground biomass and N content of winter wheat (Mistele & Schmidhalter, 2008a; Prey & Schmidhalter, 2019), the REIP index based on data from tractor-mounted sensing was used as an indicator of the spatial variability in plant biomass. In addition, there are reliable yield algorithms based on REIP (Maidl et al., 2019) that have been validated in several studies (Mittermayer et al., 2021, 2022; Schuster et al., 2022, 2023); thus, the REIP map was used to verify the rel. BMP map from satellite-sensing data.

Methods of determining soil properties

Georeferenced soil samples (SOC, TN, PCAL, KCAL, pH and texture) were collected after the grain harvest between 2018 and 2022 from the investigated. Eight soil samples at a depth of 0.3 m were taken at a maximum radius of 0.5 m around a georeferenced sampling point and combined into a composite sample. The distribution pattern of the georeferenced points was ‘systematically random’ (Thompson, 2002). SOC and TN were analyzed with a C/N Analyzer (DIN ISO 10694, 1996), KCAL and PCAL were determined by the CAL method (VDLUFA, 2012), and pH (VDLUFA, 2016) was measured with the CaCl2 method (VDLUFA, 2016). The soil texture was determined by the feel method.

Statistical and geostatistical analysis

The geostatistical data analysis was performed using R (RStudio Version 2022.12.0). The spatial resolution of the collected data varied greatly; thus, all data had to be transferred to the same raster grid and the same resolution of 10 × 10 m. Georeferenced soil samples, combine harvester data and sensor data were interpolated using ordinary kriging (Oliver & Webster, 2015). To reduce distortions of the output due to errors in the reflectance datasets, satellite imagery, combine harvester data and sensor data, outliers were cleared from the dataset using a twofold standard deviation. Before conducting ordinary kriging, for each dataset, a semivariogram was created. A semivariogram is the variance in the data according to distance classes and shows the spatial relationship of the variable with increasing separating distance (spatial autocorrelation effect). According to the distribution of the data variance and the distance classes of the semivariogram, a model was fitted. Depending on the fitted model, the data are weighted in the kriging neighborhood to predict values between sample points (Oliver & Webster, 2015). Ordinary kriging was performed for all datasets using the same target raster grid, which was based on the field boundaries of the fields.

Using the same target raster grid enabled the performance of a correlation analysis of all data points of the available datasets. Based on the Pearson correlation coefficients (r), the relationship between the digital variables and soil parameters was analyzed. The R libraries ‘tiff’, ‘rgdal’, ‘rgeos’, ‘gstat’, and ‘raster’ were used for spatial analyses and loading vector or raster files. The correlation coefficients were classified as very strong (r > 0.9), strong (0.9 > r > 0.7), moderate (0.7 > r > 0.5), weak (0.5 > r > 0.3), or very weak (r < 0.3) according to (Mittermayer et al., 2021).

Results

Results of field A1

Spatial variation in the yield pattern of winter wheat

On field A 1, was is investigated (a) whether the relative biomass potential of winter wheat determined via the NDVI shows similar distribution patterns in two cultivation years (2018 and 2020), (b) whether two-year relative biomass potential maps (2018 and 2020) have a higher accuracy and informative value than singleyear maps, (c) how close the relationships between the relative biomass potential and measured or digitally determined yields are and (d) how close the relationships between the relative biomass potential and soil parameters are.

The relative BMP maps of winter wheat in 2018 (Fig. 2a) and 2020 (Fig. 2b) showed similar yield patterns, although the weather in 2018 was clearly different from that in 2020, especially the temperature and rainfall distribution in spring (Online Resource 1). Accordingly, the relative BMP map of 2018 and 2020 (Fig. 2c) also showed a similar pattern of high- and low-yield zones.

Fig. 2
figure 2

Interpolated maps of the spatial distribution (in a 10 × 10 m grid) of the relative biomass potential of winter wheat (a–c), the winter wheat yield determined with the combine harvester (d) and by biomass sampling (e), the REIP determined by spectral reflectance measurements (f) and the soil parameters soil organic carbon (SOC) (g), total nitrogen (TN) (h) and sand content (i) at field A1

Full size image

The yield patterns derived with other digital methods (combine yield sensing system, Fig. 2d, tractor-mounted sensor system, Fig. 2f)) and the ground truth data (biomass samples, Fig. 2e) matched well with the pattern of the relative BMP maps.

The biomass potential according to the rel. BMP map of winter wheat (2018) ranged from 87.8% to 110.6% (Table 2). The absolute grain yield measured with the combine harvester yield sensing system (2018) averaged 8.6 t ha−1 and varied from 5.4 t ha−1 (62.8%) to 11.3 t ha−1 (131.4%). The absolute grain yield determined with biomass samples (2018) was higher than the yield measured at the combine, with an average of 9.9 t ha−1, varying from 6.3 t ha−1 (63.6%) to 12.9 t ha−1 (130.3%). In purely mathematical terms, this means that a change by one percent relative BMP corresponds to an absolute yield increase or decrease by approximately 3%.

Table 2 Descriptive statistic of the plant parameters (field A1)
Full size table

The soil parameters SOC, TN and texture (sand content) also showed considerable variability and clearly visible zones within the field (Fig. 2g–i and Table 3). While SOC and TN showed an almost equal distribution within the field, the sand fraction showed an inverse relationship. In areas with high sand content, SOC and NT contents were low, and conversely, in areas with low sand contents, SOC and NT contents were high. The relative BMP maps showed a similar distribution pattern as the soil parameters SOC and TN. Accordingly, the yield potential is higher when the SOC and TN contents are higher, while the yield potential is lower when the sand content is higher.

Table 3 Descriptive statistics of the soil parameters (field A1)
Full size table

Correlations between the plant and soil parameters

The rel. BMP of 2018 and 2020 were strongly correlated (r = 0.77), while the combination of 2018 and 2020 was very strongly correlated with the individual BMPs (2018, 2020) (r = 0.95 and 0.93) (Table 4). The yield measured by the combine harvester yield sensing system was moderately correlated with the rel. BMP (r = 0.57 to 0.61). The correlations of the biomass yield (r = 0.43 to 0.53) and the tractor-mounted reflectance measurements (r = 0.47 to 0.64) with the rel. BMP were weak to moderate. The relationship with the rel. BMP was greater in 2018 and in both years combined (2018 & 2020). The correlations between the relative BMP and the soil parameters SOC (r = 0.37 to 0.46), TN (r = 0.39 to 0.47) and AWC (r = 0.47 to 0.49) were closer than those to the other measured soil parameters; there was a negative relationship to the sand content (r = − 0.29 to 0.34).

Table 4 Correlation matrix (r) (field A1)
Full size table

Results of field B1

Spatial variation in the yield pattern of different crop types at different growth stages

In field B1, the pattern of the rel. BMPs of winter wheat (2018), winter barley (2019) and canola (2020) was compared at different growth stages (Fig. 3). The rel. BMP had a similar pattern at the different development stages of the analyzed crops. However, there were more or less clear differences between the relative BMP patterns of winter wheat, winter barley and canola. With the exception of canola at GS 70 (Fig. 3i), there were nonetheless stable areas in the field that had higher or lower biomass potential across crop types (Fig. 3a–h), even though weather conditions varied considerably from 2018 to 2020 (Online resource 2).

Fig. 3
figure 3

Interpolated maps of the spatial distribution (in a 10 × 10 m grid) of the relative biomass potential of winter wheat (a–c), winter barley (d–f) and canola (g–i) at crop-specific growth stages at field B1

Full size image

The variation in the rel. BMP (Online Resource 8) was initially greater for winter wheat at GS 32 and GS 39 (Fig. 3a and b) as well as for canola at GS50 and GS 60 (Fig. 3g and h) and then decreased at GS 65 and GS 70 (Fig. 3c and i). The variation in the relative biomass potential of winter barley was consistent throughout the different growth stages (Fig. 3d–f).

The biomass potential according to the rel. BMP of winter wheat (2018) ranged from 72.8 to 122.5% in growth stage GS 39 and from 88.5 to 108.4% in GS 65; thus, the variability is stage dependent. The absolute grain yield measured with the combine harvester yield sensing system (2018) averaged 10 t ha−1 and ranged from 6.2 t ha−1 (62.0%) to 13.6 t ha−1 (136.0%). A change in relative BMP by one percent is equivalent to a change in grain yield by approximately 1.5% (GS 39) to over 3% (GS 65).

For canola, the relative BMP at GS 50 ranged between 62 and 136%. At GS 60, the variation decreases (84% to 116%). After flowering (GS 70), the variation in the rel. BMP was even lower (95% to 106%).

The soil properties in field B1 varied greatly (Online Resource 9), e.g., SOC content from 1.2 to 2.1% DM, TN content from 0.11 to 0.23% DM, sand content from 16.9 to 35.2%, PCAL content from 3.4 to 12.4, and AWC from 15.9 to 22.0%. Thus, the study area was characterized by substantial differences in yield-relevant soil properties, particularly regarding nutrient supply and water retention capacity. In contrast, there were only slight differences in elevation (471 m and 487 m) across the area.

Correlations between the plant and soil parameters

The correlations between the rel. BMP at different growth stages of the crops were strong to very strong (Table 5). The strongest correlations were found for winter wheat (r = 0.94 to 0.97). The correlations between the growth stages of winter barley were strong (r = 0.74 to 0.90), and for canola, they were moderate to strong (r = 0.58 to 0.79). Between the different crops, the correlations ranged from r = 0.20 to 0.77 and were predominantly moderate. Correlations between the rel. BMP of winter wheat and winter barley were strongest at GS 39 (r = 0.65) and GS 65 (r = 0.68). The rel. BMP of canola correlated best at GS 50 with winter wheat at GS 39 (r = 0.64) and with winter barley at GS 65 (r = 0.77).

Table 5 Correlation matrix (r) (field B1)
Full size table

The highest correlation of the relative BMP of winter wheat with the grain yield of winter wheat measured by the combine harvester yield measurement system was found at GS 39 (r = 0.64). Surprisingly, the correlation between the winter wheat yield measured in 2018 and the rel. BMP of the winter barley and canola was in some cases closer than to the rel. BMP of winter wheat (e.g., r = 0.79 for canola, GS 60). This also indicates that the yield zones are very stable over the years.

The spectral reflectance measurements of the tractor-mounted sensor system correlated strongly with the rel. BMP of winter wheat during GS 32 and GS 39 (r = 0.86). This demonstrates that the satellite-based rel. BMP estimated based on the NDVI leads to a similar yield differentiation as the REIP determined by the tractor-mounted sensor system. The correlations of the rel. BMP with the soil parameters were weak to very weak. The highest correlations between soil parameters (SOC, TN, silt) and the rel. BMP of winter wheat were found at GS 65 (r = 0.40, r = 0.42 and r = 0.54). No relationships were found between elevation and rel. BMP.

Results of field C1

Spatial variation in the yield pattern of different crop types

The fields C1 and C2, which were used for model validation under practical conditions, had slightly less measurement data and results than the fields (A1 and B1) of the test stations. Nevertheless, in principle, the same analyses could be carried out, and correlations could be investigated. The rel. BMP maps of the individual crops from 2018 to 2022 (Fig. 4a–e) showed a similar yield pattern, although weather conditions were different in these years, including droughts, wet conditions and heavy rainfall (Online resource 3). The winter wheat patterns of 2018 and 2020 matched well with the rel. BMP pattern of canola 2019 (Fig. 4a, b and e). The patterns of winter barley and corn were similar to those of winter wheat, winter barley and corn, but with some differences (Fig. 4c and d). However, the main parts of the high- and low-yield zones coincided for all crops. The map of multiyear rel. BMP (Fig. 4f) showed the integrated results of the rel. BMP of all crops grown in the crop rotation.

Fig. 4
figure 4

Interpolated maps of the spatial distribution (in a 10 × 10 m grid) of the relative biomass potential of winter wheat 2018 (a), of canola 2019 (b), of winter barley 2020 (c), of corn 2021 (d), of winter wheat 2022 (e) and of the multiyear relative biomass potential 2018–2022 (f) and the soil parameters soil organic carbon (SOC) (g), total nitrogen (TN) (h) and pH (i) on field C1

Full size image

All rel. BMP of the crops showed a similar spatial variability, with minima of 91.8 to 94.5% and maxima of 105.2 to 106.9% (Online Resource 10). The rel. BMP of winter wheat (2022) and the yield derived from the tractor-mounted sensor system had a similar variation and ranged from 94.0 to 105.4% (rel. BMP) and from 8.9 t ha−1 (94.2%) to 10 t ha−1 (106.3%) (grain yield). Therefore, a 1% change in rel. BMP corresponds to an absolute change in yield by 1.07%.

The spatial distribution of the rel. BMP of all crops was consistent with the spatial distribution of the soil parameters SOC and TN, indicating higher yield potential in zones with higher SOC and TN contents (Fig. 4g and h). PH had a slightly different distribution pattern (Fig. 4i).

There was great variation in the soil parameters. The maps of SOC and TN had very similar distribution patterns. The values are shown as follows: SOC content 1.1–2.1% DM, TN content 0.12–0.22% DM, PCAL content 3.1–31.7 mg (100 g−1), KCAL content 7.0–25.3 mg (100 g−1) and pH 6.0–7.2 (Online Resource 11).

Correlations between the plant and soil parameters

The correlations between the rel. BMPs of the different crops were moderate to strong (Table 6). The highest correlation between years and crops was found between the rel. BMP of winter wheat 2018 and BMP of winter wheat 2022 (r = 0.86). The highest correlations between the rel. BMP of the crops were found between winter wheat and canola (r = 0.77 and 0.84), and the lowest correlations were found between corn and winter barley (r = 0.61) as well as between canola and winter barley (r = 0.60). The rel. BMP over all crops correlated strongly to very strongly with the individual rel. BMP of the crops (r = 0.80 to 0.92).

Table 6 Correlation matrix (r) (field C1)
Full size table

The correlations between the rel. BMP and the yield derived from the tractor-mounted sensor were moderate (r = 0.54 to r = 0.66). The highest correlation (r = 0.66) was found between REIP and the multiyear rel. BMP.

The rel. BMP maps correlated weakly to moderately with the soil parameters. The correlations with SOC and TN were highest with corn (r = 0.71 and 0.66), winter barley (r = 0.68 and 0.65) and over all crops (r = 0.66 and 0.62). The relationships between the rel. BMP of winter wheat (2018 & 2022) and SOC and TN were similar (r = 0.49 to 0.53). PCAL correlated moderately over all crops with the rel. BMP (r = 0.50 to 0.64). The highest correlation with P was found with the rel. BMP over all crops (2018–2022). The correlations with KCAL (r = 0.19 to 0.44) and pH (r = − 0.05 to 0.38) were very weak to weak.

Results of site C2

Spatial variation in the yield pattern of different crop types

Field C2 was highly variable in terms of rel. BMP (Fig. 5). A comparison of the distribution pattern of the yield potential of the individual crops revealed some differences (Fig. 5a–e). The rel. BMP maps of corn (2018), soybeans (2019), winter wheat (2020) and corn (2022) matched well in most parts of the field (Fig. 5a–c and e). However, the map of rel. BMP of winter barley (2021) showed a different distribution pattern (Fig. 5d). The multiyear rel. BMP indicated a stable yield pattern (Fig. 5f).

Fig. 5
figure 5

Interpolated maps of the spatial distribution (in a 10 × 10 m grid) of the relative biomass potential of winter wheat 2018 (a), of soybeans 2019 (b), of winter wheat 2020 (c), of winter barley 2021 (d), of corn 2022 (e) and of the multiyear relative biomass potential 2018–2022 (f) and the soil parameters soil organic carbon (SOC) (g), total nitrogen (TN) (h) and pH (i) on field C2

Full size image

The variation in the relative yield potential of corn (2018) and soybeans was low (6.98% and 7.73% variation) and ranged from 96.1 to 103.0% and from 95.7 to 103.5%, respectively (Online Resource 12). In contrast, the relative yield potential of winter barley (2021) and corn (2022) varied by a far greater margin (15.08% and 18.25%) and ranged from 90.9% to 106.0% and 91.3% to 109.6%, respectively. The variability in the relative BMP of winter wheat (2020) was not as high (93.5% to 105.5%) as the variability in the yield derived from the tractor-mounted sensor system (90.3% to 108.8%). An approximately one percent change in relative BMP corresponds to a 1.5% change in absolute yield.

The soil parameters SOC, TN and pH also showed considerable heterogeneity and distinct zones within the field (Fig. 5g–i). SOC and TN were almost evenly distributed, while pH had a different distribution pattern. The comparison of the rel. BMP maps, the multiyear rel. BMP maps and soil maps (SOC and TN) in particular, indicated that the relative yield potential is higher in areas where the SOC and TN contents are higher. The soil parameters varied considerably: SOC content 1.4–1.8% DM, TN content 0.13–0.19% DM, PCAL content 2.8–5.9 mg (100 g−1), KCAL content 3.2–14.1 mg (100 g−1) and pH 6.0–6.9 (Online Resource 13).

Correlations between the plant and soil parameters

The correlations between the crop-specific rel. BMP maps ranged from very weak to moderate (r = 0.10 to 0.53) (Table 7). Winter wheat (2022) and winter barley correlated almost equally to the other crops (r = 0.40 to 0.53). However, there was no match between corn (2022) and winter barley (2021) (r = 0.10). Corn (2018) and corn (2022) were also weakly correlated (r = 0.21). The rel. BMP map of corn (2022) only matched well with winter wheat (2020). The multiyear rel. BMP map correlated moderately to strongly with the individual crop-specific maps, while soybeans (2019) and winter wheat (2020) had the highest correlations with the multiyear map (r = 0.75 and 0.80).

Table 7 Correlation matrix (r) (field C2)
Full size table

The yield derived from the tractor-mounted sensor correlated strongly with the rel. BMP of winter wheat (2020) (r = 0.73) and the multiyear rel. BMP (r = 0.81) and correlated moderately with the rel. BMP of soybeans (2019) and corn (2022).

The rel. BMP maps, with the exception of winter barley (2021), were well correlated with the soil parameters SOC and TN. The correlations ranged from r = 0.50 to 0.73. The highest agreement with the maps of SOC and TN was observed for the multiyear rel. BMP. PCAL, KCAL and pH were only weakly correlated with rel. BMP maps.

Discussion

Discussion of methods

This study describes a new method for satellite-based remote sensing of plant-specific biomass yield patterns to determine yield zones, involving multiple observation dates, crop-specific evaluation and consideration of relevant development stages. The aim was to achieve a high level of accuracy with a simple, cost-effective and, in principle, large-scale application. The relative biomass potential is an indicator for multiyear stable and homogeneous yield zones. Areas with a higher relative biomass potential indicate a higher yield potential within the field, and information about the absolute yield level is not given by this method. For absolute yield estimation, more input data and more complicated models are needed (Klompenburg et al., 2020). However, for many applications in agriculture, knowledge of relative yield differences is sufficient, e.g., georeferenced soil sampling, measures to improve soil properties in areas with low rel. BMP, variable-rate applications of organic fertilizers, variable-rate seeding, variable-rate irrigation, and detection of unproductive areas (Bökle et al., 2023; Schuster et al., 2023).

The relative BMP maps serve as a suitable alternative to established methods based on sensor measurements or direct yield measurements for the delineation of yield zones. For instance, tractor-mounted multispectral sensor systems have sophisticated technology and methodology as well as advanced algorithms that have been sufficiently validated (Mittermayer et al., 2021, 2022; Schuster et al., 2022, 2023; Stettmer et al., 2022b); nonetheless, they are expensive and therefore not widely used in practice. They require expert knowledge and are not suitable for large-scale applications.

Yield mapping with a yield sensing system on the combine harvester is the most common method of delineating yield patterns and yield zones in modern agriculture (Chung et al., 2016). Many modern combine harvesters are equipped with yield sensing systems; however, these have a high potential for error, especially if calibration is omitted (Bachmaier, 2010; Fulton et al., 2018; Kharel et al., 2019). Yield recording with combine harvesters can provide accurate values; however, they cannot be relied upon, especially with data of unknown origin (Mittermayer et al., 2021, 2022; Stettmer et al., 2022b).

Yield measurement using a plot combine harvester and georeferenced biomass sampling is not applicable in agricultural practice and is only suitable for scientific applications. In order to derive high-quality yield maps, a sufficient number of samples must be obtained. Sample distribution as well as sample density may affect the performance of the spatial interpolation (Li & Heap, 2011); thus, the quality of the interpolated map may be insufficient if the number of samples is too low. Since the biomass samples are manually cut, they are subject to errors due to plot selection or cutting area (Petersen et al., 2005). Therefore, the biomass samples are also not the "true" values. Previous studies have shown that yields determined with biomass samples were almost always higher than yields determined with other digital and nondigital methods (Mittermayer et al., 2021, 2022; Stettmer et al., 2022a).

To conclude, the only suitable data sources for plant-based delineation of yield zones are either sensor-based or satellite-based systems. The large-scale availability of satellite data is favorable; provided that the data quality is as high as that of sensor data and that there is no loss of information. There are also satellite-based models (Promet) and commercial applications that can be used to estimate the absolute yields of crops (e.g., winter wheat) (Bach & Mauser, 2018; Hank et al., 2015). Independent validations have shown that this method is well suited for the delineation of yield zones; however, in terms of absolute yield estimations, significant yield deviations can occur (Mittermayer et al., 2021; Schuster et al., 2023; Stettmer et al., 2022b).

Discussion of results

The results at field A1 show that the rel. BMP method provides well interpretable and reliable data. The singleyear relative BMP matched well with the yield patterns of winter wheat derived from the combine harvester and derived from ground truth data (tractor-mounted sensor system and biomass samples) (Fig. 2).

The results of site B show few variations in the rel. BMP pattern during the growth period of the individual crops (Fig. 3); however, only suitable growth stages are shown. For instance, very early growth stages with low biomass development were unsuitable, as was the flowering stage in canola, which does not allow any differentiation (not shown in the figure). Therefore, the accuracy of the rel. BMP (which is based on the vegetation index NDVI) is mainly dependent on the choice of the correct growth stages according to the crop type. If incorrect growth stages are chosen, e.g., too early stages, the spatial distribution of the yield pattern will not be represented correctly, as the reflection of bare soil will lead to disturbances due to the sensitivity of the NDVI to bare soil (Mistele, 2006; Rondeaux et al., 1996). In particular, for plants that are grown in split rows, such as corn or soybeans, a growth stage in which row closure is reached must be selected. For winter wheat and winter barley, the flowering stage is one of the most suitable growth stages (Maidl et al., 2019; Prücklmaier, 2020; Spicker, 2017). This is not the case with canola. The yellow-coloured flowers lead to distortions in the yield pattern. Thus, growth stages before flowering are more suitable for canola (Spicker, 2017). For soybeans, the growth stages V5 (fifth trifoliate, R2 (full bloom) and R5 (beginning of seed filling) were selected and performed well in the yield pattern compared to the other crops. According to Andrade et al. (2022), the NDVI derived from Sentinel 2 and Landsat 8 at growth stages V5 and R2 was promising for predicting soybean grain yield.

For the delineation of yield zones, multiyear data are crucial since weather conditions are variable from year to year. The yield level in wet years can deviate from the yield level in dry years, especially if weather extremes (drought, excess water) occur (Eck et al., 2020; Gammans et al., 2017; Sjulgård et al., 2023), meaning that the yield pattern also varies. In dry years, the effects of soil properties (e.g., available water capacity and texture) on yield are more pronounced than in wet years, particularly at sites with limited water availability (Godwin & Miller, 2003; Lawes et al., 2009; Taylor et al., 2003). The relative BMP approach is a multiyear approach, so the variability in weather conditions is accounted for (Figs. 3–5).

SOC and TN are long-term stable soil parameters (Wiesmeier et al., 2019) and one of the main reasons for soil-related yield variability (Mittermayer et al., 2021, 2022; Schuster et al., 2022, 2023). As shown in the results, the relationships between the multiyear rel. BMP maps and SOC and TN were always well correlated, indicating that by this approach, the effects of soil properties on yield are well reflected (Tables 6 and 7).

Conclusion, outlook and further research

The new methodology for determining the relative BMP described in this paper can contribute to the extended application of precision farming technologies. The method is an important step towards the delineation of yield zones. The main yield zones (low-yield and high-yield zones) are well mapped using this method. The main low-yield and high-yield zones are mainly determined by soil parameters (e.g. soil texture, available water capacity) and are therefore rather stable in the long term. Influenced by the amount of precipitation, these zones are sometimes more and sometimes less important. Crop-specific patterns can vary from year to year, as the crops have different nutrient and growth requirements and there can be year-specific differences in terms of weed and disease pressure influencing the yield pattern derived from remote sensing technologies. Nonetheless, the rel. BMP approach can help to identify areas with low yield potential to manage them accordingly. This can include an adaptation of inputs to the low yield potential (e.g., seeds, mineral and organic fertilizers) or extensification by converting arable land into permanent grassland or biotopes as a contribution to environmental protection and nature conservation (Kvítek et al., 2009; Münier et al., 2004). Conversely, the BMP approach can be used to reliably delineate high-yield zones. Corresponding to the higher yield potential, these areas have higher nutrient uptake by plant stands, which must be compensated for by appropriate fertilization to avoid a decline in soil fertility (Hartemink, 2006; Sun et al., 2020). Further applications can be developed if the rel. BMP approach is consistently improved to determine both relative and absolute yield potentials. In this study, the crops grown and the respective stages of development were known. For a general scalability of this method to larger areas and regions, an algorithm that recognizes crop types and development stages dependent on the crop types must be implemented (Goldberg et al., 2021). We will address these questions in further research work. The examination of the hypothesis led to the following results:

Hypothesis 1:

The consideration of the crop type and crop-specific growth stages increases the accuracy in deriving yield zones compared to methods that do not include this information.

The results of the application of the rel. BMP method show that the choice of correct growth stages has a decisive influence on the accuracy of the derivation in yield zones.

Hypothesis 2:

The yield patterns of the respective crops are similar if the correct crop-specific growth stages are selected.

Similar yield patterns were observed among the crops, not only between winter cereals with similar nutrient and growth requirements, but also for soybeans (field C1 and partly field C2). In field C2, the distribution pattern of corn (2022) differed from that of the other crops. The correlations between the rel. BMP of the individual crops were moderate to strong (r = 0.80–0.92) in field C1. In field C2, the correlations were very weak to moderate (r = 0.1–0.51).

Hypothesis 3:

When deriving yield zones, satellite data almost reach the same accuracy as data from tractor-mounted sensor systems.

The correlations between the rel. BMP and REIP derived from the tractor-mounted sensor were (a) strong to moderate (r = 0.56–0.86), (b) the yield patterns were similar (Fig. 2), and (c) similar correlations between rel. BMP and REIP to key soil factors SOC and TN were found (r = 0.46; 0.51 at field A1, r = 0.51; 0.76 at field C1, and r = 0.50; 0.56 at field C2).

Hypothesis 4:

The spatial distribution of crop-specific relative biomass yield potential is closely related to the spatial distribution of ground truth data (e.g., biomass yield, soil properties influencing yield).

The multiyear rel. BMP was moderately to strongly correlated with SOC (r = 0.62; 0.68) and TN (r = 0.64; 0.73) (sites C1 and C2).

Thus, Hypotheses 1, 3 and 4 are confirmed. Hypothesis 2 can only be accepted with reservations.

Further research

A new methodology was successfully tested in this study. However, further validation and optimization of the BMP algorithm under completely different soil, climate and farming conditions −, e.g., large fields (> 50 ha) with extreme heterogeneity due to soil genesis −, is needed. In addition, further crop-specific validation must be carried out at different yield levels, e.g., under the conditions of organic farming.

In addition, an algorithm for estimating absolute yields will be developed by using AI methods and possibly other indices or other satellite-based spectra.