Early Modeling of the Upcoming Landsat Next Constellation for Soybean Yield Prediction Under Varying Levels of Water Availability

Abstract

The upcoming Landsat Next will provide more frequent land surface observations at higher spatial and spectral resolutions that will greatly benefit the agricultural sector. Early modeling of the upcoming Landsat Next products for soybean yield prediction is essential for long-term satellite monitoring strategies. In this context, this article evaluates the contribution of Landsat Next’s improved spectral resolution for soybean yield prediction under varying levels of water availability. Ground-based hyperspectral data collected over five cropping seasons at the Brazilian Agricultural Research Corporation were resampled to Landsat Next spectral resolution. The spectral dataset (n = 384) was divided into calibration and external validation datasets and investigated using three strategies for soybean yield prediction: (1) using the reflectance from each spectral band; (2) using existing and new vegetation indices developed based on three general equations: Normalized Difference Vegetation Index (NDVI-like), Band Ratio Vegetation Index (RVI-like), and Band Difference Vegetation Index (DVI-like), replacing the traditional spectral bands by all possible combinations between two bands for index calculation; and (3) using a partial least squares regression (PLSR) model composed of all Landsat Next spectral bands, in comparison to PLSR models using Landsat OLI and Sentienel-2 MSI bands. The results show the distribution of the new spectral bands over the most prominent changes in leaf reflectance due to water deficit, particularly in the visible and shortwave infrared spectrum. (1) Band 18 (centered at 1610 nm) had the highest correlation with yield (R² = 0.34). (2) A new vegetation index, called Normalized Difference Shortwave Vegetation Index (NDSWVI), is proposed and calculated from bands 19 and 20 (centered at 2028 and 2108 nm). NDSWVI showed the best performance (R² = 0.37) compared to traditional existing and new vegetation indices. (3) The PLSR model gave the best results (R² = 0.65), outperforming the Landsat OLI and Sentinel-2 MSI sensors. The improved spectral resolution of Landsat Next is expected to contribute to improved crop monitoring, especially for soybean crops in Brazil, increasing the sustainability of the production systems and strengthening food security in Brazil and globally.

1. Introduction

The upcoming Landsat Next, the new generation of Landsat satellites, will provide more frequent observations of the land surface, with higher spatial resolution and a greater number of spectral bands. As a result, it is expected to make a greater contribution to global agricultural monitoring, as crop systems require detailed information to identify heterogeneous areas at key phonological stages. In addition, the increased spectral resolution of 26 spectral bands [1,2] is expected to facilitate the identification of spectral signatures of specific crop traits and stresses, therefore improving crop management practices, minimizing potential yield losses, and leading to increased sustainability of production systems.

The improved resolutions (temporal, spatial, and spectral) of Landsat Next are expected to play a critical role in soybean monitoring in Brazil. Despite being the largest soybean producer in the world, Brazil´s production is constantly threatened by adverse weather conditions, especially drought [3]. In the 2023–2024 cropping season, 45.7 million hectares were dedicated to soybean production, which represents an increase of 3.8% compared to the previous cropping season [4]. However, crop production and yield were reduced due to either drought or excessive rainfall, which illustrates the spatial variability of weather conditions faced in the main production areas [4]. Specifically, in the 2023–2024 cropping season, the state of Paraná, one of the largest soybean producers in Brazil, had a soybean yield 17.9% lower than in the 2022–2023 cropping season due to low levels of water availability to support crop development [4]. In Brazil, drought periods impair approximately 30% of soybean production [5], which translates into a financial loss of USD 79 billion over 38 years, reaching over USD 14 billion during the 2021–2022 cropping season, when severe drought condition was observed in southern Brazil [6,7].

For decades, soybean yields have been estimated through remote sensing data at satellite-based [8,9,10], unmanned aerial vehicle (UAV)-based [11,12,13], and field-based levels [14,15,16], highlighting its potential for accurate, rapid and nondestructive assessment. At the field level, the hyperspectral data serve as a potential benchmark for spectral analysis, identifying essential absorption features related to crop traits, to be later transferred to multispectral sensors [17].

Thus, hyperspectral data have been successfully used to simulate multispectral data, particularly from Landsat and Sentinel sensors, for species discrimination [18], water monitoring [19], soil moisture retrieval [20], land cover mapping [21], crop nitrogen prediction [22] and aboveground biomass estimation [23,24]. This approach allows spectral information collected under controlled conditions, such as in experimental plots, to be later incorporated into orbital data to aid in model calibration, validation, and extrapolation. In soybean, satellite multispectral band simulation has proven effective for nitrogen retrieval, cultivar identification, and grain yield prediction [10,25,26].

Given that the global population is expected to exceed 9 billion people in 25 years [27], yield improvement is essential to mitigate the need to convert non-agricultural land to cropland and meet sustainable production standards [28]. Therefore, near real-time yield monitoring subsidizes public and corporate policies related to food security and food prices, storage and transportation infrastructure, and agricultural credit and insurance programs at the national and global levels [29]. In this scenario, satellite technology will certainly be constantly improved to provide superior performance compared with the technology currently available. Therefore, assessing the performance of the upcoming satellite technology shall assist policymakers in delineating future social and economic strategies regarding crop production.

Thus, early modeling of the upcoming Landsat Next products for soybean water status monitoring and yield prediction is essential for long-term satellite monitoring strategies, allowing model parameters to be fine-tuned prior to satellite operations and shortening the time gap between product release, data collection, and model development. Added to that, investigating the future spectral bands from Landsat Next might contribute to the design of new vegetation indices with the potential of enhancing the prediction capabilities of crop traits. Although many vegetation indices are widely used for crop monitoring, their fixed-band pairing and pre-existing formulas do not always translate the key spectral features and agronomic diagnosis of crop development [30].

Early modeling of the upcoming Landsat Next will not only help determine the full potential of improved bandwidths but will also serve the soybean production sector with accurate and affordable spatial datasets for agricultural management practices at regional scales [24]. Landsat Next simulated bands through ground-based hyperspectral data have proven efficient for crop residue characterization [31,32].

This article aims to evaluate the contribution of Landsat Next spectral bands for soybean water status monitoring and yield prediction under varying levels of water availability. The specific objectives are: (1) develop a PLSR model composed of key spectral bands from Landsat Next for soybean yield prediction; (2) investigate the correlation between Landsat Next spectral bands and soybean yield; and (3) evaluate existing and new vegetation indices from Landsat Next spectral bands for soybean yield prediction. The hypothesis investigated here is that the improved Landsat Next spectral bands could provide valuable information on crop water status with the potential to be correlated with yield values.

2. The Landsat Next Mission

The new generation of Landsat satellites, expected to be launched in late 2030, will continue the largest and longest-lasting space-based Earth land surface observation program that began 50 years ago with Landsat 1 and 2 [33]. Now, the new superspectral 26-band (Figure 1) will increase temporal, spatial, and spectral resolution to improve sustainable decision-making and fill information gaps for more accurate land surface monitoring, including water and aquatic systems monitoring, forest, crop and soil management, climate dynamics and mineral assessment [1].

The Landsat Next constellation consists of three identical satellites operating synchronously, with a 6-day revisit time, which is particularly important for agricultural applications because it will increase the probability of cloud-free imagery during key phenological stages, crop dynamics, and land use change [1,34]. The improved spatial resolution will provide spectral data in the visible, near-infrared, and shortwave infrared from 10 to 20 m and 60 m for atmospheric and thermal bands [1]. For agricultural applications, the improved spatial resolution, especially for the optical bands, will enhance the ability to monitor within-field crop variability and help delineate management zones.

The spectral bands of Landsat Next are compatible with the heritage bands of previous Landsat missions, have synergy with the Sentinel-2 satellite, and are also centered in wavelengths not yet explored by long-term satellite missions with information freely available to end users. Hence, Landsat Next is designed to study phenological signatures of vegetation for modeling crop development, biotic and abiotic stresses, and grain yield [1].

Table 1. Landsat Next spectral bands, spatial resolution, and wavelength range.

Band Number	Spectral Band	Spatial Resolution (m)	Wavelength Range (nm)
1	Violet	60	402–422
2	Coastal/Aerosol	20	433–453
3 *	Blue	10	457.5–522.5
4 *	Green	10	542.5–577.5
5 *	Yellow	20	585–615
6 *	Orange	20	610–630
7 *	Red 1	20	640–660
8 *	Red 2	10	650–680
9 *	Red Edge 1	20	697.5–712.5
10 *	Red Edge 2	20	732.5–747.5
11 *	NIR Broad	10	784.5–899.5
12	NIR 1	20	855–875
13	Water Vapor	60	935–955
14 *	Liquid Water	20	975–995
15 *	Snow/Ice1	20	1025–1045
16 *	Snow/Ice 2	20	1080–1100
17	Cirrus	60	1360–1390
18 *	SWIR 1	10	1565–1655
19 *	SWIR 2a	20	2025.5–2050.5
20 *	SWIR 2b	20	2088–2128
21 *	SWIR 2c	20	2191–2231
22	TIR 1	60	8175–8425
23	TIR 2	60	8425–8775
24	TIR 3	60	8925–9275
25	TIR 4	60	11,025–11,575
26	TIR 5	60	11,725–12,275

* Bands investigated in the present research.

details the spectral bands, spatial resolution, and wavelength range of the Landsat Next sensor.

3. Materials and Methods

3.1. Experimental Area

The field experiments (Figure 2) were conducted at the experimental farm of the National Soybean Research Center (Embrapa Soja), a branch of the Brazilian Agricultural Research Corporation, located in the municipality of Londrina, Paraná State, Southern Brazil (23°11′37″S and 51°11′03″W. 630 m above sea level), in the 2016–2017, 2017–2018, 2018–2019, 2022–2023 and 2023–2024 cropping seasons. All plots received the same crop management practices, including fertilizing, following the soybean production technologies [35], and the only source of variation within the experiment was the water availability to the plants.

The experimental area has an Udox Oxisol [36] with a water-holding capacity of 75 mm and a climate classification according to Köppen as Cfa climate (i.e., subtropical climate), with an average temperature in the hottest month higher than 22 °C and rainfall concentrated in the summer months, without a defined dry season [37,38], although water deficit is often observed in the summer months, coinciding with the soybean cropping season.

Following a split-plot model in a randomized complete block design with four blocks, four water treatments were distributed in the field plots (Figure 2): irrigated (IRR, receiving rainfall and, if necessary, irrigation with a soil water matric potential between 0.03 and 0.05 MPa); non-irrigated (NIRR, receiving only rainfall); water deficit induced in the vegetative stage (WDV); water deficit induced in the reproductive stage (WDR).

Table 2. Sowing dates and periods of inducement of water deficit at the vegetative and reproductive periods (expressed in days after sowing—DAS) during 2016/2017, 2017/2018, 2018/2019, 2022/2023, and 2023/2024 cropping seasons.

Cropping Season	Sowing	Water Deficit Induced at Vegetative Stages	Water Deficit Induced at Reproductive Stages	Harvesting
2016–2017	19 October 2016	From 37 DAS to 54 DAS	From 54 DAS to the harvesting period	From 116 DAS
2017–2018	18 October 2017	From 33 DAS to 62 DAS	From 62 DAS to the harvesting period	From 139 DAS
2018–2019	16 October 2018	From 41 DAS to 64 DAS	From 64 DAS to 90 DAS	From 129 DAS
2022–2023	25 October 2022	From 41 DAS to 76 DAS	From 76 DAS to 106 DAS	From 125 DAS
2023–2024	19 October 2023	From 30 DAS to 54 DAS	From 54 DAS to 83 DAS	From 109 DAS

shows the sowing dates and the periods of water deficit induced in the vegetative and reproductive stages.

In each cropping season, five soybean genotypes were distributed among the subplots. A total of 15 genotypes with drought tolerance genes and different responses to water deficit were evaluated in the five cropping seasons: 1Ea15, 2Ha11, 2Ia4, BR16, and BRS 184 in 2016–2017 and 2017–2018; Ea2939, 3Ma2, BRT18–0089, BRS 283 and BRT18–0201 in 2018–2019; M5947, BRS1064, BRS1061, M6410 and BRS 539 in 2022–2023 and 2023–2024, contributing to increase the genetic variability of the data collected.

Rainout shelters prevented WDV and WDR plots from receiving rainfall during the vegetative and reproductive phases, contributing to an increase in gravimetric soil moisture variability across treatments. During the period of water deficit induction in one treatment (WDV or WDR), the other treatment remained watered by rainfall. Shelters automatically covered the plots when the weather station located within the experimental area registered rainfall higher than 0.1 mm and uncovered plots when rainfalls ceased. Vertical concrete barriers (buried up to a depth of 90 cm) installed around the perimeter of the plots prevented lateral movement of the water from the outside into the plot soil. The irrigation schedule implemented in the IRR treatment during the 2016–2017, 2018–2019, and 2023–2024 cropping seasons (

Table 3. Irrigation schedule in the 2016/2017, 2018/2019, and 2023/2024 cropping seasons.

2016–2017 Cropping Season			2018–2019 Cropping Season			2023–2024 Cropping Season
DAS	Quantity (mm)	Duration (minutes)	DAS	Quantity (mm)	Duration (minutes)	DAS	Quantity (mm)	Duration (minutes)
24	14.4	60	52	14.4	60	22	12.0	50
29	4.8	20	53	14.4	60	47	10.8	45
30	7.2	30	57	11.5	48	48	6.0	25
31	9.6	40	58	5.7	24	55	12.0	50
34	4.8	20	59	5.7	24	56	7.2	30
35	4.8	20	61	8.4	35	60	14.4	60
36	4.8	20	66	2.9	12	61	8.4	35
37	4.8	20	106	7.2	30	63	9.6	40
38	14.4	60	109	8.4	35	64	10.8	45
			114	11.5	48	69	4.8	20
			115	2.9	12	70	4.8	20
			116	8.4	35	71	4.8	20
			119	4.8	20	75	4.8	20
						76	4.8	20
						77	4.8	20
						78	4.8	20
						81	9.6	40
						88	4.8	20
						89	4.8	20
						90	4.8	20
						91	4.8	20
						105	7.2	30
						106	8.4	35
						109	7.2	30
						111	7.2	30

) indicates 69.6, 106.2, and 183.6 mm of irrigation, respectively, which contributed to an increase in gravimetric soil moisture variability among treatments.

3.2. Weather Monitoring, Soil Moisture and Yield

The induction of water deficit and irrigation (when applied) resulted in different levels of water availability on the experimental plots, indicated by the climatic water balance according to Thornthwaite and Mather [39] for each experimental treatment of each cropping season (Figure 3). The climatic water balance was calculated from weather data (air temperature, relative air humidity, and precipitation) monitored by the weather station located within the experimental area. Additional information on weather monitoring during the evaluated cropping seasons can be found in Sibaldelli and Farias [40,41,42] and Sibaldelli et al. [43,44].

From Figure 3, it is possible to observe the efficiency of the irrigated treatment (IRR) in neutralizing the natural water deficit in the five cropping seasons monitored. In addition, the severe natural water deficit in the 2018–2019 and 2023–2024 cropping seasons is highlighted. Also, the DHV and DHR treatments proved to be efficient in promoting or increasing the water deficit during the period when the plants were deprived of rainfall (Table 2).

Figure 4 shows the effect of experimental treatments and levels of water availability on soil moisture (0–20 and 20–40 cm depths), monitored in all plots by gravimetric analysis at two growth stages: at the transition from vegetative to reproductive stage, and the R5 phenological stage towards the maturity. Growth stages were assessed weekly from emergence to maturity, according to Fehr and Caviness [45]. After fulfilling the assumptions of the analysis of variance (ANOVA), soil moisture data were subjected to ANOVA, and means were compared by the Tukey test (p ≤ 0.05) using Sisvar software version 5.8 [46].

Soil moisture assessed at the transition from vegetative to reproductive stage (Figure 4a) showed lower values for WDV, except in the 2017–2018 cropping season (0–20 cm), when the favorable natural weather conditions, as shown in Figure 3, did not affect soil moisture in this treatment. For the irrigated treatment, an increase in soil moisture was observed during the 2018–2019 (0–20 and 20–40 cm) and 2023–2024 (20–40 cm) cropping seasons. A greater influence of the different weather conditions (Figure 3) on soil moisture was observed in the R5 phenological stage (Figure 4b) for WDR. In the 2023–2024 cropping season, the severe natural water deficit condition (Figure 3) also affected soil moisture in the NIRR plot and prevented WDV from recovering its water availability as observed in the other cropping seasons.

Grain yield was calculated and corrected for 13% grain moisture, according to Equation (1):

$$\mathrm{G}\mathrm{Y}={\displaystyle \frac{\left(100\mathrm{H}\mathrm{G}\mathrm{M}\right)}{\left(100\mathrm{D}\mathrm{G}\mathrm{M}\right)}}\times \mathrm{H}\mathrm{G}\mathrm{W}\times {\displaystyle \frac{\mathrm{10,000}}{\mathrm{H}\mathrm{P}\mathrm{A}}}$$

(1)

in which GY is the grain yield (kg ha⁻¹); HGM, the harvested grain moisture (%); DGM, the desired grain moisture (%); HGW, the harvested grain weight (kg); and HPA, the harvested plot area (m²). Harvested grain moisture was measured using the G810 grain moisture meter (Gehaka Inc.—São Paulo, Brazil).

Differences in soil moisture (Figure 4), influenced by differences in water availability (Figure 3), resulted in variability in soybean yield, as shown in Figure 5. Upon meeting the assumptions of the analysis of variance (ANOVA), the yield data were submitted to ANOVA, and the means were compared by the Tukey test (p ≤ 0.05) using Sisvar software (Version 5.8).

Figure 5 and its association with Figure 4 and Figure 3 demonstrate the role of the experimental treatments in promoting differential conditions for plant growth. In 2017–2018 and 2022–2023, the nonessential irrigation resulted in similar conditions for IRR and NIRR treatments. The lack of statistical differences between IRR and NIRR in 2016–2017 is related to the low amount of irrigation performed at the initial phases of crop development (Table 2) and the similar weather conditions between them (Figure 3), leading to similar soil moisture (Figure 4a). Soybean yield from WDR plots proved to be severely impaired by water deficit. In turn, WDV plots demonstrated recovery from the water deficit induced during the vegetative stages, reaching a similar yield to IRR and NIRR plots in the 2016–2017 and 2017–2018 cropping seasons and a similar yield to NIRR in the 2018–2019 cropping season. The absence of yield differences in plants under water deficit during the vegetative stages (and rain watered during reproductive stages) compared to plants under good water availability was observed by Nogueira and Nagai [47], Rolla et al. [48], Carvalho et al. [49] and Crusiol et al. [50]. In 2022–2023, the onset of the rainout shelters in WDR, and consequently rain watering of WDV, occurred when plants were at the R3 phenological stage, undermining the flowering and pod development of WDV plants. Therefore, WDV and WDR produced similar yields. In 2023–2024, after the onset of the rainout shelters in WDR, and consequently rain watering of WDV, a severe water deficit was observed for rain-fed treatments (NIRR, WDV, and WDR—Figure 3), impairing the recovery of WDV plants from the water deficit induced.

3.3. Spectral Data Acquisition and Processing

Leaf hyperspectral reflectance was collected in the central leaflet of the fullest expanded third trifoliate leaf from the top using the FieldSpec 3 Jr spectroradiometer (Analytical Spectral Devices, Boulder, CO, USA). The ASD hyperspectral sensor has a spectral resolution of 3 nm between 350 and 1400 nm and 30 nm resolution between 1400 and 2500 nm (Figure 6). The output spectra are given in single bands of 1 nm intervals and 2151 contiguous spectral bands. Each spectral reading was the result of the average 20 internal automatic spectral readings, minimizing noise in the obtained data.

Leaves were assessed using the plant probe device (Figure 6c), connected to the FieldSpec by a one-meter bare fiber, preventing illumination interferences of adjacent targets and ensuring pure leaf reflectance spectra collection without noises, scattering, or attenuations from the interaction between electromagnetic energy and atmospheric water vapor and not requiring the use of spectral filters for noise removal and data smoothing [51]. An internal 99% reflectance board (Spectralon ^®, New Hampshire, US) was used as the reflectance standard, and a 1% reflectance blackboard was used as the opaque one. Spectral data from four plants was collected in each sub-plot and then averaged, resulting in the data used for spectral analysis.

A total of 1536 leaf reflectance samples were collected at the R5 phenological stage [45], resulting in 384 spectral samples used for data analysis, as described in

Table 4. Description of days of spectral assessment (expressed in days after sowing—DAS) and the number of spectral samples in the 2016–2017, 2017–2018, 2018–2019, 2022–2023, and 2023–2024 cropping seasons.

	Cropping Season
	2016–2017	2017–2018	2018–2019	2022–2023	2023–2024
DAS	89	96	94	92	82
Spectral samples	64	80	80	80	80

. The R5 phenological stage represents the most suitable timing for yield prediction in soybean crops under variant water status conditions through remote sensing data [8,10,15,16,29,52]

3.4. Statistical Analysis

Following the methodological flowchart displayed in Figure 7, the obtained hyperspectral dataset was resampled to the spectral resolution at the full width at half maximum of Landsat Next spectral bands at 10 and 20 m spatial resolution, excluding band 2, designed for water quality assessments. The visible, near-infrared and shortwave infrared bands investigated were: 3, 4, 5, 6, 7, 8, 9, 10, 11, 14, 15, 16, 18, 19, 20 and 21 (Table 1). Additionally, to ascertain the full potential of Landsat Next concerning currently available spectral datasets, we also resampled the hyperspectral data to the Sentinel-2 (bands 3, 4, 5, 6, 7, 8, 11, and 12—[53] and Landsat 8 (and 9) satellites (bands 2, 3, 4, 5, 6 and 7—[54]).

The spectral data were submitted to principal component analysis (PCA) to explore the qualitative spectral differences among soybean water status, represented by the experimental treatments in each cropping season. To prospect the quantitative potential of Landsat Next for yield prediction, we evaluated the possibility of generating a multiyear yield prediction model. For that purpose, spectral data from 2016 to 2017, 2017 to 2018, 2018 to 2019, 2022 to 2023, and 2023 to 2024 and the correspondent yield data were pooled into the same dataset, totaling 384 samples.

Spectral data (n = 384) were randomly split into two subsets: calibration/cross-validation (containing 75% of data—288 spectral samples) and external validation (containing the remaining 25% of the data—96 spectral samples), used to test the developed models. Figure 8 presents the yield (Figure 8a) and the spectral (Figure 8b) data distribution, demonstrating the equality between calibration/cross-validation and external validation datasets.

We investigated the correlation between each Landsat Next spectral band and the yield and also evaluated the capability of existing and new vegetation indices for yield monitoring. Additionally, Landsat Next, Sentinel-2, and Landsat 8 (and 9) datasets were also submitted to partial least squares regression (PLSR) to ascertain the development of a soybean yield prediction model composed of reflectance derived from all Landsat Next spectral bands. PCA and PLSR were performed by The Unscrambler^® (CAMO Software—Norway).

3.4.1. Principal Component Analysis (PCA)

Principal component analysis was applied to the spectral dataset from each cropping season to investigate how effectively the spectral response of soybean plants under variant water status could be explained and clustered. PCA is a data mining method that reduces the number of variables to be analyzed (spectral bands) by transforming them, using a covariance matrix composed of all spectral bands, into a new group of variables named principal component (PC). PC is the linear combination of all analyzed spectral bands to explain their variance within the dataset, and its score represents the percentage of explained variance. The first principal component (PC1) carries the most spectral information on data variance; the second principal component (PC2) carries the residual information of PC1; the third principal component (PC3) carries the residual information of PC1 and PC2 (PC1 + PC2), and so on [55,56].

3.4.2. Landsat Next Spectral Bands

In the calibration/cross-validation dataset (containing 75% of data—288 spectral samples), reflectance from each of the 16 Landsat Next spectral bands was correlated to soybean yield. The coefficient of determination (R²) from the linear regression between yield and the correspondent reflectance was used to select the outstanding spectral band, which was then applied to the external validation dataset (containing 25% of data—96 spectral samples).

The fit quality of the linear model using the outstanding spectral band was assessed by the coefficient of determination (R²—Equation (2)) and the root mean square error (RMSE—Equation (3)) between observed and predicted values. Although RMSE provides a robust metric for fit quality assessment of a specific model, it does not consider data variability and distribution and, therefore, the comparison among models developed from different datasets (calibration and external validation) is weakened. To overcome this limitation, the RMSE is often normalized by the standard deviation, representing, thus, the Residual Prediction Deviation (RPD—Equation (4)) [57]. However, RPD might be affected by extreme values. Therefore, we analyzed the ratio between the interquartile distance and RMSE, named Ratio of Performance to Interquartile Distance (RPIQ—Equation (5)), to assess the fit quality of the adjusted models, considering both the root mean square error and data distribution [58]. The higher the RPIQ, the better the model fit.

$${\mathrm{R}}^2=1-{\displaystyle \frac{\mathrm{S}\mathrm{S}\mathrm{R}}{\mathrm{S}\mathrm{S}\mathrm{T}}}=1-{\displaystyle \frac{\sum {({\mathrm{y}}_{\mathrm{i}}-{\widehat{\mathrm{y}}}_{\mathrm{i}})}^2}{\sum {({\mathrm{y}}_{\mathrm{i}}-\overline{{\widehat{\mathrm{y}}}_{\mathrm{i}}})}^2}}$$

(2)

where $SSR$ is the Sum Squared Regression and $SST$ is the Sum of Squares Total.

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{1}{\mathrm{n}}}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{X}}_{\mathrm{i}}-{\mathrm{x}}_0\right)}^2}$$

(3)

where $n$ is the number of variables, and $X_i$ is the observed value and $x_0$ the predicted value.

$$\mathrm{R}\mathrm{P}\mathrm{D}={\displaystyle \frac{\mathrm{S}\mathrm{D}}{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}}$$

(4)

where $SD$ is the standard deviation.

$$\mathrm{R}\mathrm{P}\mathrm{I}\mathrm{Q}={\displaystyle \frac{(\mathrm{Q}3-\mathrm{Q}1)}{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}}$$

(5)

where Q3 is the third quartile, and Q1 is the first quartile.

3.4.3. Landsat Next Vegetation Indices (VIs)

We assessed all possible combinations between two spectral bands to calculate vegetation indices. Many of the calculated VIs already exist and are currently used to monitor crop traits [59]. However, the procedure adopted here adjusts the bandwidth to the spectral resolution of Landsat Next and explores the capability of the improved spectral bands of Landsat Next to generate new vegetation indices to monitor soybean yield.

Vegetation indices were calculated based on three general formulas: Normalized Difference Vegetation Index (NDVI-like—Equation (6)), Band Ratio Vegetation Index (RVI-like—Equation (7)), and Band Difference Vegetation Index (DVI-like—Equation (8)), as per:

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}={\displaystyle \frac{{\mathrm{B}\mathrm{a}\mathrm{n}\mathrm{d}}_{\mathrm{x}}-{\mathrm{B}\mathrm{a}\mathrm{n}\mathrm{d}}_{\mathrm{y}}}{{\mathrm{B}\mathrm{a}\mathrm{n}\mathrm{d}}_{\mathrm{x}}+{\mathrm{B}\mathrm{a}\mathrm{n}\mathrm{d}}_{\mathrm{y}}}}$$

(6)

$$\mathrm{RVI}={\displaystyle \frac{{\mathrm{B}\mathrm{a}\mathrm{n}\mathrm{d}}_{\mathrm{x}}}{{\mathrm{B}\mathrm{a}\mathrm{n}\mathrm{d}}_{\mathrm{y}}}}$$

(7)

$$\mathrm{DVI}={\mathrm{B}\mathrm{a}\mathrm{n}\mathrm{d}}_{\mathrm{x}}-{\mathrm{B}\mathrm{a}\mathrm{n}\mathrm{d}}_{\mathrm{y}}$$

(8)

where Band_x and Band_y can be replaced by bands 3, 4, 5, 6, 7, 8, 9, 10, 11, 14, 15, 16, 18, 19, 20 and 21. For the calculated VI, the position of the two spectral bands for soybean yield prediction in Equations (6)–(8) will determine the tendency (positive or negative) for the correlation between yield and the VI. If the position of the two spectral bands is changed, the R² remains the same, but a different tendency (positive or negative) will be observed.

A total of 360 vegetation indices was calculated—120 for each general formula. In the calibration/cross-validation dataset (containing 75% of data—288 spectral samples), vegetation indices were correlated with soybean yield. The coefficient of determination (R²) from the linear regression between yield and the correspondent VI was used to select the outstanding vegetation index from each general formula, which was then applied to the external validation dataset (containing 25% of data—96 spectral samples), and evaluated by the R² between the predicted and observed values of yield and the RPIQ.

3.4.4. Partial Least Squares Regression (PLSR)

Partial Least Squares Regression (PLSR) is a multivariate regression model used to correlate spectral data to a trait variable (e.g., yield) and is especially advantageous to cope with the multi-collinearity usually found in spectral data [60] when the number of the predictor variable (e.g., spectral bands) is larger than the number of response variables (e.g., yield). PLSR generates a new dataset of orthogonal base vectors (PLSR factors) that account for most of the variation in a trait variable, delivering a linear model consisting of waveband scaling coefficients to transform full-spectrum data [61].

Up to date, many machine-learning algorithms are being used for crop traits estimations through remote sensing. However, no single machine-learning algorithm fits all data, highlighting the need to choose specific algorithms for each prediction scenario [62]. Added to that, most machine-learning algorithms used for prediction tasks must be fine-tuned, which is not always carried out, undermining their full potential [57]. Hence, we chose the PLSR model based on its well-established methodology for model adjustment and based on our previous research for soybean yield prediction using spectral data [10,15]. The comparison among different machine-learning models for yield prediction is beyond the goals of the present research.

PLSR models were developed (based on the reflectance from the 16 spectral bands from Landsat Next, from the six bands from Landsat-8 and the eight bands from Sentinel-2) considering the optimal number of orthogonal base vectors (a key process that deeply affects the model prediction capacity—[51], represented by the lowest value of root mean square error (RMSE), highest coefficient of determination (R²), and value of Bias close to zero through the “leave-one-out” cross-validation method [63]. Before the model development, reflectance spectra were normalized by subtracting the mean reflectance from the actual reflectance at each wavelength, enabling the comparison among the waveband scaling coefficients. Outliers and homogeneity of the spectral data assessment were assessed by the Leverage and Hotelling’s T2 tests.

In a second stage, aiming at assessing the prediction performance within each cropping season, PLSR models were also developed individually for each cropping season using reflectance correspondent to the Landsat Next spectral bands and evaluated at the calibration and cross-validation steps. Additionally, aiming at assessing the influence of genetic variability on the prediction performance, PLSR models were developed for each soybean cultivar within each cropping season.

Fit quality of the PLSR models was assessed by the coefficient of determination (R²) between observed and predicted values and the Ratio of Performance to Interquartile Distance (RPIQ), considering both the root mean square error and data distribution and enabling the comparison among PLSR models developed individually for each cultivar, cropping season and the multiyear model.

4. Results and Discussion

4.1. Soybean Spectral Behavior as Affected by Water Availability and Contribution of Landsat Next Spectral Bands for Water Status Monitoring

Figure 9 shows the spectral response of soybean plants under IRR and WDR treatments during the 2023–2024 cropping season. Plants under water deficit demonstrated higher reflectance across the Vis-NIR-SWIR spectrum than plants under adequate water availability (Figure 9a), which is consistent with Damm et al. [64]. It should be noted that reflectance differences across wavelengths vary significantly due to the biochemical properties and structure of leaves [65,66].

The primary effect of crop water status on soybean leaf reflectance is demonstrated across SWIR wavelengths and is directly related to the absorption of radiation by water [67,68]. Thus, SWIR reflectance has been shown to have a strong negative relationship with leaf water content and, consequently, crop water status [69,70,71]. The secondary effects are demonstrated across Vis and NIR wavelengths and cannot be explained by water absorption of radiation alone. At visible wavelengths, leaf pigments and the absorption of photosynthetically active radiation play a remarkable role in light absorption, which is directly impaired under water deficit. Consequently, lower gas exchange rates occur, which consequently impairs CO₂ influx. Under these conditions, reflectance increases at visible wavelengths to avoid degradation of the photosynthetic apparatus due to excess light [67,72,73,74,75,76]. NIR wavelengths are affected by light scattering along the mesophyll, which depends on internal leaf structures such as cell wall width, intercellular air spaces, and the amount of mesophyll per unit of leaf area within the mesophyll [66,67,77,78].

Figure 9b presents the percentage of reflectance increasing from WDR treatment compared to IRR, with the delimitation of the correspondent Landsat Next spectral bands (color bars). The lowest increments in reflectance were observed across NIR wavelengths, demonstrating larger increments across Vis and SWIR wavelengths. The Landsat Next spectral bands demonstrate even distribution across the visible spectrum, covering the peaks of reflectance increasing observed at 580 nm with improved spectral resolution (bands yellow, orange, and red 1 bands) and at 702 nm (red edge band 1). Across the four near-infrared bands, a very small percentage of increasing reflectance was observed. The increasing reflectance peaks across SWIR wavelengths at 1400 and 1990 nm are located on the atmospheric absorption window, and no spectral bands are centered in those spectral intervals. Landsat Next bands 18, 19, 20, and 21 seem to be located on the most conspicuous changes in soybean reflectance across the SWIR spectrum. Specifically, bands 19, 20, and 21 are in a spectral interval represented by a single band in Landsat OLI and Sentinel-2 MSI sensors, which might significantly contribute to soybean water status monitoring.

Aiming at exploring the contribution of improved spectral resolution from Landsat Next for soybean water status monitoring, Figure 10 presents the principal component analysis for each cropping season.

In the 2016–2017 cropping season (Figure 10a), WDV and WDR treatments demonstrated different distributions between them and from IRR and NIRR in the scatterplot of PC1 and PC2, highlighting the effects of water deprivation in their spectral behavior. In 2017–2018 (Figure 10b), under the absence of irrigation and natural water deficit, only WDR could be consistently clustered. Similar behavior was observed in the 2018–2019 cropping season (Figure 10c) when WDR was consistently clustered. Among the other three treatments, NIRR presented an intermediate distribution between WDV and IRR. In the 2022–2023 cropping season (Figure 10d), the absence of natural water deficit prevented the clustering of treatments, with a discrete trend of clustering for WDV plants subjected to water deficit during vegetative stages, flowering, and pod development stages. In the 2023–2024 cropping season (Figure 10e), the severe natural water deficit led to similar spectral behavior between NIRR and WDV, with a trend of clustering out IRR treatment. Only WDR was consistently clustered. The PCA results evidence the potential of Landsat Next for clustering soybean water status, especially in plants that experienced water deficit periods, such as WDV and WDR treatments.

4.2. Spectral Modeling for Soybean Yield Prediction Based on Landsat Next Reflectance

Based on the differential spectral behavior in soybean plans under variant levels of water availability, this section aims to explore the quantitative potential of Landsat Next for yield prediction.

4.2.1. Landsat Next Spectral Bands

Figure 11 displays the coefficient of determination from the linear regression between reflectance from each Landsat Next spectral band and soybean yield in the calibration dataset (n = 288). The lowest coefficients of determination were observed across the near-infrared spectrum, while the highest R² was observed across SWIR wavelengths. Among the visible spectrum, band 9 presented the highest R².

The lowest R² values for the NIR spectral bands are consistent with the percentage of increasing reflectance shown in Figure 9b, where NIR bands demonstrated little increase in reflectance due to water deficit. Similarly, Landsat band 9 showed an increasing reflectance peak (702 nm, Figure 9b) and the highest R² among the visible bands. While NIR bands have low sensitivity to changes in water status, band 9, centered at 705 nm, may be related to photosynthetic activities that are directly affected by water deficit. Photosystem I, in which the reaction center preferentially absorbs light at 700 nm, plays a crucial role in electron transport during photosynthesis, thus interfering with the process of photosynthetic activity inhibition.

Despite competitive values of R² across SWIR bands, bands 18 and 19 presented the highest values, reaching R² up to 0.36 and 0.35, respectively. The results are relevant since the highest correlation between soybean yield and reflectance was observed in new spectral bands from Landsat Next, such as bands 19, 20, and 21, not available in the Landsat OLI sensor. Albeit bands 9 and 18 are available in Sentinel-2 MSI sensor, band 19 will be exclusively available in Landsat Next.

When applying the linear regression between yield and reflectance from band 18 (outstanding spectral band in Figure 11) to external samples (external validation dataset, n = 96), a positive correlation was achieved between the observed and predicted yield values (Figure 12), with, however, low R² (0.34). It should be emphasized the similar data distribution in the calibration and external validation datasets, demonstrating that both yield data (Figure 8a) and spectral data (Figure 8b) have similar distribution in both datasets, strengthening the calibration and prediction steps. Considering the need for exploring the full potential of Landsat Next improved spectral resolution, processing, and analytical methods, such as vegetation indices, might increase the prediction ability.

4.2.2. Landsat Next Vegetation Indices (VIs)

We assessed all possible combinations between two spectral bands to calculate vegetation indices based on Normalized Difference Vegetation Index (NDVI-like—Equation (6)), Band Ratio Vegetation Index (RVI-like—Equation (7)), and Band Difference Vegetation Index (DVI-like—Equation (8)), which comprises the majority of existing vegetation indices and expands the horizons for optimized new vegetation indices based on the improved spectral resolution from Landsat Next for yield prediction.

Figure 13 presents the coefficient of determination from the linear regression between soybean yield and all possible combinations between two spectral bands to calculate vegetation indices based on normalized difference vegetation index (NDVI-like—a), ratio vegetation index (RVI-like—b), and difference vegetation index (DVI-like—c) in the calibration dataset (n = 288).

A similar pattern was observed for NDVI (Figure 13a) and RVI (Figure 13b), with low R² when using Vis-Vis, Vis-SWIR, and NIR-NIR bands as input data. Intermediate values of R² were observed for Vis-NIR band combinations. It should be noted that the input using band 9 combined with bands 11, 15, and 16, which reinforces its contribution to yield monitoring, as observed in Figure 9b and Figure 11.

Outstanding band combinations were observed when using SWIR bands as input data. The highest R² (0.45) was observed combining bands 19 and 20 as input for index calculation. This result is higher compared to vegetation indices traditionally used for yield monitoring, such as the traditional NDVI (using bands NIR and red bands 8 and 11 from Landsat Next, Rouse et al. [79]), NDII1 (using bands 11 and 18 from Landsat Next, Hardisky et al. [80]), NDII2 (using bands 11 and 21 from Landsat Next, Hardisky et al. [80]) and the NDRE (using bands 9 and 11, Gitelson and Merzlyak, [81]), which delivered R² equal to 0.17, 0.17, 0.23 and 0.30, respectively.

Although competitive results were obtained for the RVI using bands 19 and 20, the normalized difference of bands 19 and 20, named here as the Normalized Difference Shortwave Vegetation Index (NDSWVI), enables data analysis through time-based on data normalization from −1 to 1.

For DVI (Figure 13c), a different pattern from NDVI and RVI was observed, with the highest coefficient of determination found when combining the red edge 1 band (band 9) with visible bands as input for index calculation. The outstanding difference vegetation index was calculated using bands 4 (green) and red-edge1 (band 9), with R² equal to 0.41.

When the linear regression between yield and the outstanding vegetation indices was applied to external samples (external validation dataset, n = 96), a positive correlation between observed and predicted yield values was obtained (Figure 14). Using NDVI and RVI indices calculated with bands 19 and 20, an improved predictive ability was observed compared to a single spectral band (Figure 12), with R² equal to 0.37 for both indices, although a higher RPIQ was observed for RVI. However, the DVI index calculated with bands 4 and 9 underperformed the use of a single spectral band (R² = 0.32) (Figure 12).

The optimized vegetation index using the improved spectral resolution from Landsat Next proved the best performance when calculated with new spectral bands. Specifically, for the NDVI calculated with bands 19 and 20, both bands are located within a spectral interval represented up to date by a single band, both in Landsat OLI and Sentinel-2 MSI sensors, highlighting the Landsat Next contribution for yield prediction.

4.2.3. Partial Least Squares Regression (PLSR)

Figure 15a presents the correlation between observed and predicted soybean yield from the calibration stage of PLSR (n = 288), developed to ascertain the potential of a soybean yield prediction model composed of reflectance derived from all Landsat Next spectral bands.

The coefficient of determination from the calibration stage of PLSR (0.71) outperformed the results obtained for a single spectral band (0.36) and the NDVI calculated with bands 19 and 20 (0.45). PLSR also proved more efficient (R² up to 0.65) than spectral bands and vegetation indices during the external validation (Figure 11 and Figure 13).

Figure 16 presents the contributions of Landsat Next spectral bands for the PLSR model, demonstrating the highest importance of bands 19 and 20, the same bands used for calculating the outstanding optimized vegetation indices (Figure 12 and Figure 13). Obtained results reinforce, once again, the contribution of spectral bands from the SWIR spectrum from Landsat Next for soybean monitoring, especially bands 19 and 20, represented within the same spectral band on Landsat OLI and Sentinel-2 MSI sensors.

4.3. Comparison of Prediction Performance Using Spectral Bands, Vegetation Indices and PLSR Models

Figure 17 gathers the performance, expressed by the R² and RPIQ, in the external validation stage of the spectral models developed using the outstanding spectral band (band 18), outstanding vegetation indices, and PLSR models using all spectral bands from Landsat Next, Landsat OLI, and Sentinel-2 MSI.

Among the strategies used for soybean yield prediction using Landsat Next spectral bands, the spectral model composed of reflectance from all spectral bands delivered superior results compared to single spectral bands and vegetation indices. Spectral bands pooled into multivariate or machine-learning regression models are recognized to present higher prediction accuracy compared to single spectral bands or vegetation indices for soybean yield [9,10,82].

The improved spectral resolution of Landsat Next was shown to outperform Landsat OLI (R² = 0.35) and Sentinel-2 MSI (R² = 0.50), reinforcing the contribution of the new spectral bands (both in amount, width, and center) of Landsat Next. The underperformance of Landsat OLI compared to the Sentinel-2 MSI sensor was also reported by Sibanda et al. [24] in predicting the aboveground biomass of grass species. Prey and Schmidhalter [17] performed all possible band combinations to calculate vegetation indices for winter wheat yield prediction and obtained better results for the Sentinel-2 sensor than the Landsat OLI.

A slight decrease in R² and RPIQ was observed when comparing Landsat Next spectral bands to hyperspectral data (data not shown), with R² equal to 0.69 and RPIQ equal to 2.79 at the external validation step. According to Inoue et al. [59], Ullah et al. [83], Yendrek et al. [61], and Crusiol et al. [84], the simultaneous use of hundreds of spectral bands from hyperspectral data has leverage over broadband multispectral data in the prediction of crop traits due to its narrow key spectral features. The Landsat Next showed competitive results between both spectral resolutions proving the contribution of its bands for yield prediction.

4.4. Perspectives for Improving Landsat Next Spectral Models for Soybean Yield Prediction

The results of PLSR are promising since soybean yield values ranging from 287 to 5398 kg ha-1 could be predicted under different water availability within each cropping season and different climatic conditions among the five cropping seasons evaluated (Figure 3). The spectral complexity of different cropping seasons could be challenging to a multiyear spectral model development for soybean yield prediction. Thus, we developed PLSR models individually for each of the five cropping seasons evaluated aiming at reinforcing the prediction performance of the Landsat Next spectral bands. The distribution of yield data within each cropping season is presented in Figure 18, revealing different patterns for cropping seasons without natural drought (2017–2018 and 2022–2023) and cropping seasons under the influence of natural water deficit (2016–2017, 2018–2019 and 2023–2024).

PLSR models (Figure 19) presented R² ranging from 0.67 to 0.92 at the calibration step and from 0.59 to 0.88 at the cross-validation step through the leave-one-out method. The Ratio of Performance to Interquartile Distance (RPIQ) assessed at the calibration step reveals that only the 2017–2018 and 2023–2024 cropping seasons (Figure 19b and Figure 19e, respectively) outperformed the multiyear PLSR model (Figure 15a), with RPIQ equals to 3.23 and 3.22 compared to 2.91 from the multiyear model. Obtained results highlight the decay of accuracy with increasing spectral complexity and strengthen the potential of spectral models composed of reflectance from all spectral bands to predict yield in cropping seasons with severe natural drought conditions, such as 2023–2024.

Additionally, we investigated the influence of genetic variability within the spectral dataset (15 different genotypes) on PLSR models. Figure 20 presents the statistics of the cross-validation step of PLSR models developed individually for each soybean genotype within each cropping season.

High variability in the accuracy among soybean genotypes within each cropping season was observed, with R², ranging from 0.19 (cultivar 2Ia4 on 2016–2017 cropping season) to 0.91 (genotype 3Ma2) on 2018–2019 cropping season, while RPIQ ranges from 0.98 (genotypes 1Ea15 on 2016–2017 cropping season) to 4.51 (genotypes M5947 on 2023–2024 cropping season).

The genetic variability is known for having a direct impact on soybean spectral response [13,85,86,87]. Thus, genotypes assessed under equal water availability might have differential physiological responses, delivering differential grain yield [48,71,88]. Consequently, differential spectral behavior is expected among soybean genetic materials, which poses a challenge to yield modeling through spectral data, as highlighted by Christenson et al. [14].

In this context, the characterization of most soybean genotypes adopted in the leading soybean production regions ensures higher stability of the spectral models. In addition, there is a need to overcome limitations in the transferability of spectral models [10,89,90], which requires enlarging the spectral dataset to gather spectral data collected over a variety of edaphoclimatic conditions for soybean production in Brazil and under different management practices, strengthening the prediction ability of Landsat Next dataset for yield prediction.

The closely related spectral behavior existing between ground- and satellite-based multispectral data, albeit differences in reflectance intensity [10,91,92,93], indicates the outstanding performance of Landsat Next improved spectral resolution and reinforces the contribution of early modeling of the forthcoming Landsat Next products for soybean water status monitoring and yield prediction for long-term satellite monitoring strategies, enabling fine-tuning of model’s parameters before the satellite operation and shortening the time gap between product release, data collection and model development.

5. Conclusions

This article evaluates the contribution of Landsat Next spectral bands for soybean water status monitoring and yield prediction under varying levels of water availability, showing the distribution of the new spectral bands over the most conspicuous changes in leaf reflectance due to water deficit, especially over the visible and shortwave infrared spectrum.

Near-infrared bands showed a low correlation with soybean yield, while shortwave infrared bands showed a better relationship with yield values, with the highest correlation for band 18, centered at 1610 nm. Despite a lower correlation, band 9 (centered at 705 nm) showed the best performance among the visible bands. New vegetation indices developed from the improved spectral resolution of Landsat Next outperformed existing traditional vegetation indices. The outstanding vegetation index, called the Normalized Difference Shortwave Vegetation Index (NDSWVI), was calculated from the normalized difference of bands 19 and 20 centered at 2028 nm and 2108 nm, respectively. Among the strategies for soybean yield prediction using Landsat Next reflectance, the PLSR model, using reflectance from all spectral bands, proved superior to the use of single spectral bands or vegetation indices (both existing and new indices), outperforming the Landsat OLI and Sentinel-2 MSI sensors.

The improved spectral resolution of Landsat Next will help improve crop monitoring, especially for soybean crops in Brazil; increase the sustainability of production systems in agronomic, environmental, economic, and social terms, as demonstrated by better management and higher yields; prevent conversion of non-agricultural land to cropland; provide greater profitability and subsidize public and corporate policies for the commercial market, agricultural credit programs, and crop insurance contracts; and strengthen food security in Brazil and globally.

Future research should address the limitations that might result from the methodology adopted in the present article, such as the impact of spatial and temporal resolution rather than spectral resolution. How the improvements in the spatial and temporal resolution will benefit the soybean yield estimation is yet to be reported. Added to that, this article focused on the assessment of yield as dependent on the plant’s water status, considering drought occurrence as one of the key factors impairing soybean yields. Therefore, the methodology described here might be expanded for characterizing pests, diseases, and nutrient deficiencies that might impact soybean yield and spectral response.