Located in the Changping District of Beijing ( Figure 1 ), Xiaotangshan National Precision Agriculture Demonstration Base (40.10°N, 116.26°E; altitude 39 m) has a typical climate of the northern winter wheat zone, with an average of 2506.5 hours of sunshine a year, an average annual temperature of 13.3°C, and an average annual rainfall of 563.8 mm. The base is used for high-quality agricultural research area relevant to large irrigated areas and high winter wheat yields.

The ability of the four common vegetation indices to predict the winter wheat yield was evaluated using four methods (LSTM, RF, GBDT, and SVR); the results of these predictions are visualized in Figure 3 . A comparison showed that, of the four models, the LSTM made the best predictions, followed by the RF model. The spatial distribution of the yield predicted by the four models roughly matched the true yield distribution: in each case, the yield was high in the middle of all the plots and low near the edges, a pattern that may have been due to human activities and the presence of trees around the plots. Overall, it was shown that the LSTM, RF, GBDT, and SVR models could be used to make estimates of the winter wheat yield that also reflected the spatial distribution.

All three metrics (R², MAE, and RMSE) showed that the LSTM model produced the best yield estimates (see Figure 4 ), with a value of 0.93 for R². The neural network architecture of the deep learning model also performed well, giving an R² of 0.886. The values of R² for the GBDT and SVR models were 0.839 and 0.573, respectively. The two tree-based models – RF and GBDT – were able to explain the yield change at least 10% better than the SVR. In contrast, although the SVR technique could effectively solve multiple collinearity problems among independent variables, it only simulates the limited relationship between input variables/features and modeling targets (i.e., grain yield), and is unable to map highly non-linear and complex relationship between variables. As reported in many previous works, deep learning methods are generally considered to be superior when the number of training samples is sufficiently large (Maimaitijiang et al., 2020; Khaki et al., 2021). This is likely due to the fact that DL often exceeds popular machine learning methods when dealing with larger sample size, complex, nonlinear and redundant datasets (LeCun et al., 2015). Our experiments verified that the network containing two LSTM layers could capture more than 90% of the yield information from the input features. Previous research has also demonstrated that LSTM model performed best through several machine learning models in winter wheat prediction. Xie and Huang (2021) demonstrated that the accuracy of the LSTM model was significantly higher than that of the 1-D CNN model due to the better ability of the LSTM model to treat time-series satellite data. However, when the amount of data is limited, the RF model has the advantages of being insensitive to outliers, nonlinearity, serial autocorrelation, and high dimensionality. For example, Cao et al. (2021) found that the performance of RF was not always worse than DL at both the county and field levels. What’s more, although LSTM yielded superior performance over RF methods, the improvement in grain yield prediction accuracy was not substantial (see Figure 4 ), one reason could be the little difference in measured gridded yield data. Future work will examine the ability of more advanced deep learning architectures (e.g., LSTM and its variants) at county scale to extract better information for winter wheat yield prediction.

The LSTM method was then used to produce estimates of the winter wheat yield based on 30-m ZY-1 02D data, 30-m Sentinel-2 data, and 10-m Sentinel -2 data. The modeling with the ZY-1 02D data used the vegetation indices NDSI, SSI, and RSI designed by ourselves as features, whereas three conventional vegetation indices – NDWI, SR, and EVI – were used with the Sentinel-2 data. The time-series consisting of 20210324, 0330, 0428, and 0501 were input into the LSTM model for training, and predictions of the yield were obtained for the entire study area (see Figure 5 ). It can be seen that, for all three types of data, the yield distribution was correctly modeled and that 10-m Sentinel-2 data best reflect the actual distribution of the yield within the study area. The spatial resolution of the results based on the other two datasets is low; as a result, the corresponding yield distribution maps are coarse and do not reflect the differences in yield between adjacent grid cells.

Overall, the estimates based on the 10-m Sentinel-2 data were found to be the most accurate, followed by those based on the 30-m ZY-1 02D data; the estimates based on the 30-m Sentinel-2 were the least accurate (see Figure 6 ). By comparing the results for the datasets with different spatial and spectral resolutions, it was found that the spatial resolution had a greater impact on the winter wheat yield estimates than the spectral resolution: this can be seen from a comparison of the results for the 10-m Sentinel-2 data and the 30-m ZY-1 02D data. The 10-m Sentinel-2 data, which was the dataset with the highest spatial resolution, performed best, capturing 91% of the yield variation. The observed data consisted of gridded data with a spatial resolution of 5 m, and the satellite data with the spatial resolution that was closest to this produced the best estimates of the yield. The result was also recognized in previous studies that hyperspectral PRISMA models was lower than the multispectral Sentinel-2 models (Marshall et al., 2022). However, multispectral bands provide coarser spectral information than hypersectral bands (Yang et al., 2021). A comparison of the results based on the ZY-1 02D and Sentinel-2 data, which have the same spatial resolution, showed that the ZY-1 02D data, which has more spectral bands, performed better, indicating that the features most important to the yield still remained after the band-by-band combination and that the narrow bands could provide relevant and accurate information about the yield. Therefore, Hyperspectral (HS)-Multispectral(MS) fusion paradigm to hyperspectral data is considered to get both advantages of high spatial and spectral resolution. Here, some of the advantages and limitations of applying hyperspectral imaging to estimates of agricultural yields are demonstrated. However, although the mathematical relevance of the hyperspeactral band selection method is maximized, the computational volume is large and the physical meaning is not clear enough, resulting in low applicability (Kong et al., 2022). Further research should be targeted towards to the implementation and evaluation of more applicable band select method in hyperspectral, making the most effective use of hyperspectral band information.

In the PFI experiment that was conducted, the importance of the four features used in the modeling could be ranked as SR > NDWI > EVI > REP. The shuffled vegetation index SR produced the largest loss of 0.4783, followed by a loss of 0.4492 for the NDWI, 0.2385 for the EVI, and 0.2371 for the REP. A larger loss value indicates a greater contribution to the results. The large SR contribution is due to the high correlation between the ratio of the red band to the NIR wavelength bands and the leaf area index, which is a good measure of the crop growth (Jordan, 1969). The NDWI indicates the amount of biostructural water contained in a crop, so the large contribution made by the NDWI indicates that water has a great influence on the accumulation of organic matter in a crop (Gao, 1996). The EVI, which is a greenness vegetation index, had less effect on the results, which may be related to the instability of the blue band due to residual atmospheric effects. The poor performance of the REP vegetation index in the modeling may be due to the fact that the spatial resolution of the red-edge bands in Sentinel-2 is 20 m, which does not match the spatial resolution of the yield data. However, the vegetation indices we selected were based on anthropogenic experience about spectral information; other vegetation indices such as Green Leaf Area Index (Duchemin et al., 2008), Crop Water Stress Index (Ghaemi et al., 2016) should also be considered. We can combine empirical and statistical methods in the selection of vegetation indices.

In most previous studies, each feature was input to a model individually to retrain the model (Cai et al., 2019). This seems an intuitive approach; however, it is not appropriate if we are interested in the feature importance of the model where all features are trained together instead of one by one. Zhang et al. (2021) evaluated six typical VIs separately for their abilities to predict maize yield using the three approaches. Compared the feature importance rank in Zhang et al. (2021), the advantage of the PFI method that we used in this study is that it outputs the performance of each feature when all of the features are input to the model together, and all interactions with other features are automatically considered. We analyzed the spectral information in combination with the environmental stresses of the crop, and the ranking results reflected the most important factors for wheat growth, so as to provide a reference for practical agricultural management: in conducting winter wheat farming, we need to focus on natural conditions such as tillage density, which is related to the leaf area index, and the amount of irrigation, which is closely related to the water within the winter wheat plant. However, besides PFI, more and more methods have recently been proposed to help users interpret the predictions of complex models, such as a unified framework for interpreting predictions named SHapley Additive exPlanations (Lundberg and Lee, 2017), and further research should be targeted towards to the comparison among multiple feature importance rank methods.

From the analysis of the correlation between the different vegetation indices and the yield based on the band-by-band combination method, we determined which bands had the strongest positive and negative correlations with the three different vegetation indices (see Table 3 ). Of the 13203 combinations calculated for the NDSI, the largest correlation with the yield was for the wavelength range 516–765 nm. This lies in the visible and NIR region and was negatively correlated with the yield with a correlation coefficient of –0.7413. For the SSI, the best combination was 671–679 nm, which was negatively correlated with the yield with a correlation coefficient of –0.7559. For the RSI, the strongest correlation was for the combination 1779–2216 nm, which was positively correlated with the yield with a correlation coefficient of 0.7539. The correlation between bands combination and yield is consistent with (Marshall et al., 2022). These three customized vegetation indices were then adopted as three features for input to the hyperspectral data model. The correlation coefficients shown in Table 3 all pass the significance test of P-value< 0.001. Importantly, from the values of the correlation between the vegetation indexes and the yield (see Figure 7 ), it was found that the combinations of visible and NIR bands were negatively correlated with the yield, whereas there was a positive correlation between the shortwave infrared narrowband combinations and the yield. And can be seen from Table 3 , the difference between the most positively correlated shortwave infrared narrowband combination and the most negatively correlated visible–NIR band combination is almost negligible (between 0.5% and 5%), which indicates that both combinations can provide equally important information for yield estimates. The combinations of bands around the maximum value in Figure 7 can almost play the same role as the chosen bands (see Table 3 ) and can also be used to build the yield estimation model when the requirements are not too strict. This conclusion confirmed the results reported in the previous studies choosing the visible and NIR bands to compute VIs (Jin et al., 2017) for yield prediction to various crops. Kong et al. (2022) used the band-by-band combination method between 450-950 nm to construct new vegetation index, and analyzed the correlation of them with LAI. In our work, we expanded band range to full bands following the recommendation of Marshall et al., 2022 to get a comprehensive use of hyperspectral information.

The high degree of correlation between the visible–NIR band combinations and the yield can be explained by the correlation between the crop growth and the chlorophyll content (Acito et al., 2022), which has also been demonstrated in previous studies (Zhang et al., 2021). However, in most studies, only bands in the range 500 to 900 nm have been used and most other useful bands have been neglected (Zhang et al., 2018; Zhang et al., 2021). In this study, we also found that the shortwave infrared bands between 1000 and 2500 nm had a highly positive correlation with the yield (see Figure 7 ), suggesting that combinations of narrow shortwave infrared bands can provide equally important information to that provided by visible–NIR bands for crop yield estimation and that the use of the two types of information might achieve even better results. In future studies, experiments using the two types of bands should be performed to determine the quantitative relationship with the winter wheat yield.

The results described above (see Figures 3 , 4 ) show that the LSTM model produced significantly better estimates of the winter wheat yield than the other three models; these results are consistent with those found by Zhang et al. (2021) and Lin et al. (2020). Compared with machine learning models, deep learning with complex neural network structures has the advantage that it processes high-dimensional data that reflects the growth and development of crops (Mu et al., 2019; Wang et al., 2020). The LSTM can learn more time-dependent information (Hochreiter and Schmidhuber, 1997). The data input to the models consisted of series of monthly data that can be used to explore various types of changes in crop growth, including the yield and other related crop parameters (Haider et al., 2019; Tian et al., 2021; Wang et al., 2022). The two tree models (RF and GBDT) proved to be better than the SVR model at estimating the yield, which also confirms the results reported in a previous study (Lin et al., 2020). The LSTM neural network provides an effective tool for building new data-driven models for regional yield estimation. Neural network modeling transforms raw input variables into high-level representations through nonlinear activation and squashing functions, which weakens the traceability and interpretability of the LSTM model (Tian et al., 2021). You et al. (2017) add a DGP to Long Short Term Memory (LSTM) network, and outperforming all the competing approaches. Future work will examine the ability of more advanced deep learning architectures (e.g., LSTM and its variants) and more kinds of data (e.g., remote sensing data and climate information) to extract better multimodal information for grain yield prediction. The introduction of an attention mechanism to capture and interpret the contribution of each time node in the time-series data to the models can be considered; in combination with knowledge of the crop phenological period, this could be used to make estimates of pre-production early yields. In subsequent studies, the use of transfer learning methods to improve the scalability of the model could also be tried; this would be similar to a method of predicting winter wheat FVC using deep transfer learning (Yu et al., 2022).

Experiments were then performed in which the Sowing, cultivation and irrigation were varied. It was found that the winter wheat accumulated the most organic matter when the seeding rate was set at 225 kg/ha (0.8 times the conventional sowing rate), rotational tillage was adopted, variable amounts of fertilizer were used, and the amount of irrigation was set at 60 mm (see Figure 8 ). It was also found that a high sowing density leads to a lower yield, probably due to the intense competition between water, fertilizer, and light making the plant less biologically productive and causing lodging (Chen et al., 2022). Figure 8B shows that the choice of rotational tillage as the tillage practice can maximize soil fertility: some studies have shown that rotational tillage practices can increase the soil porosity and improve the nutrient quality, thus increasing crop yields (Nie et al., 2015). From Figure 8C it can be seen that the use of variable fertilization promotes yield improvement because it meets the nitrogen demand of winter wheat throughout the growing season and allows the crop to maintain a more reasonable canopy structure for photosynthesis even after flowering (Jiang et al., 2015). Figure 8D shows that the greater the amount of irrigation, the more water is absorbed by the crop due to osmotic pressure regulation; this increases the soil water storage and improves the drought tolerance of the wheat and the yield. The above meteorological variables are closely related to the crop growth process and directly affect the yield. Currently, as chemical fertilizer prices are rising sharply and many places are experiencing water shortages, advanced studies that will lead to the application of precise amounts of fertilizer and irrigation so that planting costs and environmental pollution can be reduced and high yields of wheat achieved are research priorities (Carberry et al., 2013). Our study quantitatively explored the effect of environmental conditions on the winter wheat yield, and the results provide data that are important to the cultivation of winter wheat in northern China.

Meanwhile, from the above experimental results, it can be seen that crop management statistics have a strong correlation with the crop yield and can be used to indicate yield changes. Management statistics models are included in crop environment models, which can be used in crop yield estimation by establishing a correlation between the crop management statistics and crop yields (Guarin and Asseng, 2022). The most commonly used statistical management model is the systemic integrated factor forecasting method developed by Chen (1992). This forecasting method predicts the annual grain yield by building a systematic model between statistical factors (irrigation, fertilizer usage, and mechanical inputs) and the crop yield. Besides management statistics models, crop environment models (Launay and Guerif, 2005) also include agrometeorological models and agronomic yield estimation models – the former use integral regression models based on meteorological factors and yields (O'Neal et al., 2002), and the latter mainly establish relationships between crop growth conditions and crop yield components, thus allowing them to predict crop yields. Predictions of wheat yields based on management statistics, meteorological data, and crop growth conditions can be highly accurate; however, these models do not apply to large areas and the values of the parameters are difficult to determine. Using satellite remote sensing, crop information can be acquired repeatedly over large areas at a relatively low cost. The combination of remote sensing data with crop growth dynamics models to predict yields has shown promise, and several studies (Cao et al., 2020; Zhang et al., 2021; Beyene et al., 2022) have shown that combining remote sensing and other factors can improve the accuracy of yield estimates. The previous studies demonstrate the tremendouspotential of remote sensing data-based crop yield prediction when employing a multimodal data fusion and deep neural network approach. Maimaitijiang et al. (2020) verified that multimodal data fusion yielded superior performance for yield prediction over single sensor data, regardless of modeling methods. Therefore, in addition to being adaptable to different remote sensing data-VIs, within-field, multi-field, and regional applications require grain yield models to cope with variation and heterogeneity in space caused by differences in soil, irrigation, fertilization and other field conditions that affect plant growth (Maimaitijiang et al., 2020). For example, Su et al (2017) integrated geographical data from the weather station in China and the SVR method to estimate crop growth at various stages. In subsequent research, we plan to collect multiple types of data, including meteorological data and crop management statistics, and combine these with remote sensing data to produced more accurate winter wheat yield estimates.

However, there remain challenges to fully understanding changes in winter wheat yields that arise from a lack of understanding of the mechanisms involved or a lack of data. If fused or integrated data with a high temporal, spatial, and spectral resolution (Cao et al., 2020) can be obtained, transferring our proposed model to a larger study area can be considered. In this study, we found that an increase in either the spectral or spatial resolution leads to an increase in the estimation accuracy ( Figure 4 ). Acito et al. (2022) deals with the problem of improving the spatial resolution of hyperspectral data from the PRISMA mission and provides a superresolved image with a spatial resolution of 10 m and the same spectral resolution as the PRISMA hyperspectral sensor. In future, following the work by Acito et al. (2022), In future, we also hope to use fused ZY-1 02D and Sentinel-2 data for yield estimation to explore how temporally, spectrally, and spatially rich data perform, in a similar way to how the fusion of ZY-1 02D and multispectral data has been used for land classification (Sun et al., 2020a). Compared with yield estimation models on county scale (Sun et al., 2020b), the yield estimation models that we developed in this study all apply at the pixel scale and may be less applicable at the larger scale of experimental fields. In subsequent studies, we will collect yield data from a large number of counties and cities, compare agricultural fields in their natural state with small experimental fields under human management, and explore the performance of the yield estimation models over spatially heterogeneous large plots. Furthermore, the approach can be tested for different crop types at different development stages and environmental conditions to evaluate the robustness.