Solar-Induced Chlorophyll Fluorescence (SIF) has emerged as a promising remote sensing index in recent years owing to its close coupling with photosynthesis, demonstrating its superior performance in agricultural monitoring and environmental assessment (Zhang et al., 2018). Compared with traditional vegetation indices such as the Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI), SIF offers more direct information about photosynthesis, making it advantageous for crop yield estimation and drought monitoring (Sun et al., 2015; Mohammed et al., 2019; Qiu et al., 2022). This study utilized a global SIF dataset from the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/) provided by Yao (Yao, 2021), which covers SIF observations on a global scale (4d, 0.05°×0.05°). In this study, the raw data were resampled using bilinear interpolation with a spatial resolution of 500 metres and the coordinate system was converted to WGS-84 to improve accuracy and applicability. To ensure the image time alignment problem. Images with the same or similar dates as the start of winter wheat sowing in EVI and LAI were selected, and the maximum value compositing (MVC) method was used to generate an 8-day time series image of winter wheat during the growing period, and cloud cover and atmospheric disturbances were effectively reduced by selecting the maximum SIF value. County-level wheat SIF values were then averaged to obtain the mean SIF for each region during the growth period, providing a more accurate assessment of photosynthesis and solid foundation for subsequent yield estimation and drought monitoring.

Two other satellite data products commonly used in crop yield estimation, including the Leaf Area Index (LAI) and Enhanced Vegetation Index (EVI), were selected. Research has indicated that LAI, representing the total leaf area per unit of land area, can be closely related to crop growth and yield accumulation, making it a key indicator of vegetation coverage. An EVI optimized to minimize atmospheric and soil noise can be widely adopted for monitoring vegetation growth and coverage. However, the original MODIS LAI (MOD15A2H, 8d, 500m) and MODIS EVI (MOD13A1-16d,500m and MOD13A1-16d,500m combined as 8d, 500m) products often contain discontinuities and noise owing to cloud and precipitation interference (https://lpdaac.usgs.gov/). To address this issue, a Savitzky-Golay (Window Size=9, Polynomial Order=2) filter (Chen et al., 2004) was applied to smooth the raw data pixel-by-pixel, effectively reducing noise while preserving primary trends. Subsequently, county-level LAI and EVI values were averaged over an 8-day time series with a spatial resolution of 500m using the winter wheat planting period as the first day to improve the accuracy and reliability of the data.

The interaction between climate change and wheat yields is complex. The climate data used in this study were derived from ERA5(1d, 0.1°×0.1°), a fifth-generation atmospheric reanalysis dataset produced by the European Center for Medium-Range Weather Forecasts (ECMWF, https://www.ecmwf.int/). ERA5 combined the model outputs with global observations to provide a comprehensive global dataset, including the 2-meter air temperature (T2m), minimum and maximum air temperatures (Tmn and Tmx), total precipitation (Pre), the 10-meter u-component wind speed (U10m), and the 10-meter v-component wind speed (V10m). For climatic data with shorter time intervals (1 day), they were converted to 8-day temporal resolution by data aggregation (8-day aggregated averaging) based on the start dates of other datasets, and the coordinate system was converted to WGS-84 to ensure consistency with other datasets. Secondly, winter wheat meteorological data were obtained by masking method using winter wheat distribution information, and finally aggregated to county-level wheat meteorological data by averaging method. The 2-meter air temperature and temperature range were crucial for understanding crop growth, as they could affect photosynthesis and growth rates. The total precipitation affected the soil moisture and water availability for wheat, and the wind components could provide insights into evaporation rates and potential extreme climate events. Analyzing these variables is essential for understanding the impact of climate change on wheat yield. Additionally, we used the Standardized Precipitation Evapotranspiration Index (SPEI) to represent drought conditions across different time zones in the study area. The SPEI data (https://spei.csic.es/index.html) provides global long-term drought information with a spatial resolution of 0.5 degrees and a monthly temporal resolution. This dataset features multi-scale properties, offering SPEI time scales from 1 month to 48 months. We utilized the SPEIbase v2.9 version, with data spanning from January 1901 to December 2022. For our study area, we obtained monthly SPEI data from 2010 to 2020 at a spatial resolution of 0.5°×0.5°. Drought classification according to SPEI is provided in Table 1 .

Reliable winter wheat yield data are essential for enhancing the accuracy of model yield estimates. This study used county-level winter wheat yield data from the National Bureau of Statistics of China (https://www.stats.gov.cn/), focusing on 10 consecutive years and selecting counties with yield levels above 3,000 kg/ha for analysis. After screening, 93 counties were included in each year. Statistical yearbook data indicated that the annual variation in the winter wheat planting area was less than 1%, with winter wheat being the primary crop from October to early June. The cropland data from the land use classification were applied as a mask for winter wheat, with the MODIS land classification product (MCD12Q1) providing yearly images and selecting land cover types based on the IGBP global vegetation classification scheme. These data were processed using the Google Earth Engine (GEE), and annual winter wheat planting distribution maps for each county were generated in combination with county administrative boundary vector files.

The experimental approach in this study was divided into three stages ( Figure 2 ): data collection and processing, model construction comparison, and yield prediction. Initially, multisource remote sensing data (including MOD15A2H, MOD12A1, MYD12A2, and MCD12Q1), climate data (EAR-5), and the novel remote sensing index SIF for Henan Province from 2011 to 2020 were collected. During data preprocessing, all data were resampled to an 8-day interval with a 500 m resolution, and the Savitzky-Golay filter was adopted to smooth the time series data(EVI, LAI). The mean of each feature was then calculated, the target area wheat data were extracted through mask processing, and all features were normalized to a range of [0, 1] to prepare them for the model input.

In the model construction and yield prediction stage, various models, including traditional algorithms such as RF and XGBoost, as well as certain deep learning models such as BO-CNN-BiLSTM (BCBL) and baseline LSTM、Transformer, were trained using the preprocessed data. Bayesian Optimization was employed to fine-tune the hyperparameters of the models and enhance their estimation accuracy. The study then performed a spatiotemporal analysis of winter wheat yield from 2015 to 2020 and assessed the predictive ability of the models, with a focus on the accuracy of predictions made 25 d prior to harvest. The prediction accuracy of each model was compared using performance metrics, such as R², RMSE, and MRE. Additionally, the role of SIF data in enhancing the accuracy of traditional remote sensing data was analyzed, with a particular emphasis on its performance under extreme climate conditions.

Bidirectional Long Short-Term Memory (BiLSTM) is an advanced version of the Long Short-Term Memory (LSTM) network designed to effectively handle time-series data while preserving memory. Based on LSTM’s solution to the vanishing and exploding gradient problems inherent in Recurrent Neural Networks (RNNs), BiLSTM can enhance time-series modeling by running two LSTM networks simultaneously. One network can process the data forward in time, and the other processes it backward, enabling the consideration of both past and future information.

The key components of LSTM include cell states and several gating mechanisms: forget, input, and output gates. These mechanisms can regulate the flow of information between different gates using a sigmoid function.

For a given time step t and input X t , each LSTM unit is defined as follows ( Figure 3 ):

(1) Forget gate:

(2) Input gate:

(3) Updating unit status:

(4) Unit state:

(5) Output gate:

(6) Hide status updates:

In this context, σ represents the sigmoid function.

BiLSTM computes the bidirectional hidden states by running LSTM layers in both the forward and backward directions along the time axis. The forward LSTM processed the sequence from the first to the last element, whereas the backward LSTM processed it in reverse order. The resulting hidden states from both directions were concatenated to form the final bidirectional hidden state.

Bayesian optimization is an effective global optimization method that is particularly suited for costly black-box objective functions. Hyperparameter optimization for machine learning models guides each iteration by incrementally constructing a probabilistic model, typically a Gaussian process, of the objective function, thereby efficiently identifying the optimal solution.

The key steps involved training an agent model f ( x ) using existing data points and applying Gaussian process regression to estimate the distribution of the objective function. Subsequently, the next sample point x n + 1 was selected for evaluation using an Acquisition Function that balances exploration (discovering new regions) and exploitation (leveraging known good regions). The acquisition function we use is Expected Improvement (EI). The Bayesian optimiser estimates the mean and uncertainty of the objective function based on the agent model (Gaussian process) and then selects the location where the Expected Improvement is the largest to be sampled, with the goal of finding the point that maximises the predictive uncertainty of the agent model, usually expressed as:

where α ( x | D n ) represents the acquisition function of the current dataset D n . By incorporating new data points in each iteration, Bayesian optimization gradually converges to the global optimal solution. Its advantage is the effective optimization of complex functions with a limited number of evaluations, making it particularly suitable for tasks such as hyperparameter tuning.

The convergence criteria for Bayesian optimization primarily rely on two aspects: optimization iterations and improvement threshold. Optimization iterations involve setting a maximum iteration count, n_iter, at which point the optimization process will terminate once the specified count is reached. In our case, we chose n_iter = 50, indicating 50 optimization iterations. Improvement threshold entails selecting the next sampling point based on the expected improvement criterion from the surrogate model. After a certain number of iterations, if the expected improvement becomes minimal, indicating limited room for further enhancement in the current region, it is considered as convergence, and the optimization process gradually ceases.

The steps for Bayesian optimization in selecting hyperparameters can be described as follows ( Figure 4 ):

Step 1: Initialization: In the initial phase, a certain number of hyperparameter combinations are randomly selected for evaluation (we set init_points = 10, choosing 10 points at random for evaluation).

Step 2: Modeling: A Gaussian process or other surrogate model is used to model the objective function, predicting its shape based on existing evaluation results.

Step 3: Sampling Point Selection: The next sampling point is chosen based on the surrogate model’s predictions and the expected improvement criterion.

Step 4: Evaluation and Update: The selected hyperparameter combination is used to evaluate the objective function (in this context, training and validating the model), and the surrogate model is updated based on the evaluation results.

Step 5: Iteration: Steps 2, 3, and 4 are repeated until the stopping criteria are met (such as maximum iterations or convergence standards).

Compared to traditional grid and random search, Bayesian optimization is more computationally efficient and converges faster. First, it reduces the number of evaluations by using the surrogate model to efficiently explore the parameter space, avoiding redundant calculations in previously evaluated regions. This typically allows it to find optimal hyperparameters with fewer evaluations. Second, it does not rely on exhaustive evaluations, as traditional methods like grid search do, which require evaluations across the entire parameter space. Instead, Bayesian optimization leverages an efficient sampling strategy for global optimization with minimal computational cost. In summary, Bayesian optimization excels at optimizing complex functions within limited evaluations, making it suitable for tasks like hyperparameter tuning. The hyperparameter optimization range for the model in this study can be found in Appendix S1 in Supplementary Data Sheet S1 .

In this study, a BO-CNN-BiLSTM (BCBL) deep learning model was developed for feature extraction from remote sensing, climate, and yield time-series data. The model consisted of an input layer, a convolutional layer (1DCNN), a pooling layer, a reshaping layer, a bidirectional LSTM layer, two fully connected layers, and an output layer and was implemented using the TensorFlow framework, as shown in Figure 5 . The input features included EVI+LAI (TS), SIF, and climate data, with new features generated through various combinations of these inputs. The target vector was the county-level yield data from 2012 to 2020. Initially, the feature combinations were normalized, and the dataset was randomly split into 80% training and 20% testing sets. Given the significant impact of hyperparameters on model performance, Bayesian global optimization was employed to optimize the hyperparameters, such as the number of units in the convolutional layer, LSTM layer, fully connected layer, and learning rate. Model testing was then conducted using optimal hyperparameter configurations.

Two machine learning models were compared, and deep learning model (LSTM、Transformer) was selected as the baseline to evaluate the accuracy of the models. Random Forest (RF) generates final predictions by constructing multiple independent decision trees and averaging their results, demonstrating strong nonlinear modeling capabilities and being widely used in bagging methods. XGBoost, as a gradient boosting-based algorithm, progressively optimized the model by iteratively constructing decision trees that fit the residuals of previous trees, achieving notable performance in regression and classification tasks. As a baseline, an LSTM model with one LSTM layer and two dense layers was chosen, while a Transformer model with MultiHeadAttention, residual connections, layer normalization, fully connected layers, and dropout was also selected. All models underwent hyperparameter tuning using Bayesian optimization, with the RF and XGBoost results averaged to determine the optimal model after five cross-validations.

To evaluate the estimation performance of the models, three metrics were adopted: coefficient of determination (R²), root mean square error (RMSE), and mean relative error (MRE). The formulas for these metrics are as follows.

where x i is the statistical yield of wheat; y i is the estimated yield of winter wheat; y ¯ i is the mean statistical yield of wheat; and n is the number of samples in the test set.

According to the model and dataset combinations outlined in Methods (2.3), SIF outperformed all other models in the single-feature datasets (including TS, climate, and SIF), surpassing both TS and climate. This superior performance was due to the close association of SIF with plant photosynthesis and its rich yield information, which is crucial for accurate crop yield estimation. In contrast, using TS and weather data alone produced similar results but with lower accuracy, likely because of their coarser resolution and mixed image elements. Additionally, comparing the climate data with the TS+climate data revealed only minimal improvements in accuracy, with both the machine learning and deep learning models demonstrating low final accuracies. This indicated that the optimization effect of climate data on TS information was limited, with the contribution of climate data being gradually absorbed by TS data as crop fertility progressed, resulting in minimal improvement. Among all the data combinations, the TS+climate+SIF combination yielded the best performance in the deep learning models, achieving a coefficient of determination (R²) of 0.81 for the BCBL model. Conversely, the TS+SIF combination delivered the best results in the machine learning models, with a coefficient of 0.73 for the XGBoost model ( Figure 6 ).

Deep learning can generally excel with large data samples owing to its multilevel architecture, which can learn the deep connections between features, unlike machine learning performing less effectively in this regard. In terms of model accuracy, XGBoost and RF demonstrated superior performance with the TS+SIF dataset. However, adding climate data slightly reduced the accuracy of these models. This decrease may be attributed to two factors: first, the climate data did not significantly enhance the TS dataset in machine learning; second, the vast amount of data and its complex interrelationships limited the ability of machine learning to conduct deep data mining, resulting in poorer performance with the TS+climate+SIF combination.

Throughout the study, the accuracy of all four models, whether machine learning or deep learning, was improved to varying extents by incorporating SIF into the TS and TS+climate data combinations ( Figure 7 ). Specifically, in the TS+climate+SIF combination, SIF enhanced the R² values of the BCBL, Transformer, LSTM, XGBoost, and RF algorithms by 0.23,0.16, 0.21, 0.23, and 0.23, respectively, while reducing the RMSE by 305.13,227.29, 246.41, 258.40, and 230.70 kg/ha. In the TS+SIF combination, SIF increased the R² values of BCBL, Transformer, LSTM, XGBoost, and RF by 0.23,0.41, 0.27, 0.37, and 0.35, respectively, and reduced the RMSE by 270.11,449.15,304.31, 398.37, and 368.40 kg/ha. These results indicated that SIF significantly complemented the yield estimation as a remote sensing index and enhanced model performance across different datasets. Consequently, BCBL and LSTM exhibiting superior performance were selected for further study.

Figure 8 presents a scatterplot comparing the estimated yields versus the official statistical yields for the BCBL and baseline LSTM models using the same dataset throughout the study period. In the TS+Climate+SIF dataset, the single LSTM model achieved an R² of 0.75, an RMSE of 715.86 kg/ha, and an MRE of 8.26% ( Figure 8B ). In contrast, the BCBL model benefiting from CNN feature extraction and bidirectional LSTM two-way temporal learning achieved an R² of 0.81, an RMSE of 616.99 kg/ha, and an MRE of 7.14% ( Figure 8A ). The BCBL model thus improved R² by 0.06, and reduced RMSE and MRE by 98.87 kg/ha and 1.11%, respectively, compared to the LSTM model. These improvements indicate that the BCBL model significantly enhances yield prediction accuracy and robustness under different data conditions.

The scatterplot shows that the yield estimates of both models are evenly distributed around the 1:1 line, displaying a clear linear relationship, which demonstrates the model’s ability to capture the trend between yield and prediction effectively in most cases. Compared with the single-layer LSTM model, the BCBL model’s scatter points are more concentrated around the 1:1 line, indicating that its predictions are closer to the actual values. Particularly in low- and high-yield areas, the BCBL model significantly reduces overestimation in low-yield regions and underestimation in high-yield regions, resulting in more accurate and stable predictions. By integrating CNN spatial feature extraction with bidirectional LSTM two-way temporal learning, the BCBL model more comprehensively captures the spatiotemporal features of the crop growth process, effectively reducing prediction errors. Further analysis of residuals shows that the absolute values of most residuals are within 750 kg/ha, with only a few counties exceeding this range, indicating the BCBL model’s high accuracy in yield estimation across most areas. However, for counties where residuals exceed 750 kg/ha, data quality issues or model underfitting may be contributing factors. Specifically, high residuals in these counties may result from insufficient inclusion of local meteorological, soil, or other environmental factors in the model, or the model’s weaker generalization ability in specific areas. These outliers indicate that although the overall performance of the BCBL model is better than that of LSTM, there is still room for further optimization in the future to better handle the regional variability problem.

In summary, the BCBL model demonstrates a clear advantage over the LSTM model in yield prediction, with smaller errors, and most of the estimated value points are near the 1:1 line, thus reflecting the trend of yield more reliably, and having an improved effect on the underestimation of high yield and overestimation of low yield.