The Ability of Sun-Induced Chlorophyll Fluorescence From OCO-2 and MODIS-EVI to Monitor Spatial Variations of Soybean and Maize Yields in the Midwestern USA

: Satellite sun-induced chlorophyll fluorescence (SIF) has emerged as a promising tool for monitoring growing conditions and productivity of vegetation. However, it still remains unclear the ability of satellite SIF data to predict crop yields at the regional scale, comparing to widely used satellite vegetation index (VI), such as the Enhanced Vegetation Index (EVI) from the Moderate Resolution Imaging Spectroradiometer (MODIS). Additionally, few attempts have been made to verify if SIF products from the new Orbiting Carbon Observatory-2 (OCO-2) satellite could be applied for regional corn and soybean yield estimates. With the deep neural networks (DNN) approach, this study investigated the ability of OCO-2 SIF, MODIS EVI, and climate data to estimate county-level corn and soybean yields in the U.S. Corn Belt. Monthly mean and maximum SIF and MODIS EVI during the peak growing season showed similar correlations with corn and soybean yields. The DNNs with SIF as predictors were able to estimate corn and soybean yields well but performed poorer than MODIS EVI and climate variables-based DNNs. The performance of SIF and MODIS EVI-based DNNs varied with the areal dominance of crops while that of climate-based DNNs exhibited less spatial variability. SIF data could provide useful supplementary information to MODIS EVI and climatic variables for improving estimates of crop yields. MODIS EVI and climate predictors (e.g., VPD and temperature) during the peak growing season (from June to August) played important roles in predicting yields of corn and soybean in the Midwestern 12 states in the U.S. The results highlighted the benefit of combining data from both satellite and climate sources in crop yield estimation. Additionally, this study showed the potential of adding SIF in crop yield prediction despite the small improvement of model performances, which might result from the limitation of current available SIF products. The framework of this study could be applied to different regions and other types of crops to employ deep learning for crop yield forecasting by combining different types of remote sensing data (such as OCO-2 SIF and MODIS EVI) and climate data.

Many studies have proved that satellite SIF products are able to track vegetation production (e.g., gross primary production, GPP) for various ecosystems [49,[59][60][61][62][63]. The integration of satellite SIF and process-based models could improve the mechanistic understanding and monitoring of crop productivity [48,64]. However, the application of satellite SIF data for estimating crop yield using the statistical approach has less been practiced, mainly owing to the low spatial resolution of the abovementioned SIF data. For example, the footprint sizes of GOME-2 and SCIAMACHY SIF are 40 km 80 km and 30 km 60 km, respectively. Within a pixel of such large sizes, several crops types might coexist, inducing a big challenge in estimating crop yield using coarse resolution SIF data in heterogeneous regions since relationships between SIF and vegetation production change with crop types [48,65,66]. Recently, Cai et al. [24] demonstrated the applicability of GOME-2 and SCIAMACHY SIF data in conjunction with climate data in predicting wheat yield in Australia, in which wheat is dominantly distributed.
Fortunately, the recent advent of the Orbiting Carbon Observatory-2 (OCO-2) [55] and TROPOspheric Monitoring Instrument (TROPOMI) [67] enables relative high-resolution SIF retrievals with footprint areas of 1.3km 2.25 km and 3.5km 7 km, respectively. At such spatial resolutions, it is possible to develop SIF-based statistical yield prediction models with the differentiation of different types of crops, such as C3 and C4 crops. However, OCO-2 SIF data has limitations for estimating crop yields, such as discrete footprints and low revisiting coverage. Therefore, it is worth investigating whether OCO-2 SIF data is applicable for monitoring crop yield at regional scales.
In this study, the ability of OCO-2 SIF to monitor the spatial variations of soybean and maize yields in the Midwest 12 states of the U.S.A. was investigated. The DNN method was employed to construct prediction models of yields with OCO-2 SIF, MODIS EVI, climate variables, and their different combinations as inputs. The performance of different models was evaluated using census data of county-level soybean and maize yields. The importance of various variables during different periods in predicting crop yields was assessed. The specific objectives of this study are: (1) to assess the ability of OCO-2 SIF to monitor the spatial variation of maize and soybean yields relative to EVI; (2) to compare the ability of remote sensing and climate data to monitor the spatial variations of corn and soybean yields in the Midwestern US; (3) to identify key periods when remote sensing and climate data are important for predicting crop yields.

Materials and Methods
The framework of this study is shown in Figure 1. Firstly, we processed the OCO-2 SIF and MODIS EVI data to generate county-level means of SIF and EVI for corn and soybean, respectively (section 2.2.3 and 2.3.1). Secondly, climate data from the Daily Surface Weather and Climatological Summaries (DAYMET) databases [68], including temperature, precipitation, and vapor pressure (VP), were processed for each county using the Google Earth Engine (GEE) platform. Lastly, DNN models were trained with census yields and above processed data (section 2.3.2).

Study Area
The study area (Figure 2) covers the 12 states of Midwestern US, including Illinois, Indiana, Iowa, Kansas, Michigan, Minnesota, Missouri, Nebraska, North Dakota, Ohio, South Dakota, and Wisconsin. This region, known as the Corn Belt, made up more than 80% of the harvest area and over 82% of the production of either corn or soybean in the US from 2015 to 2017 ( Figure S1). The growing season in this region generally spans from April to September. More detailed information on the soil, climate, and practices of management in this region can be found in [69].

Satellite Data
OCO-2 SIF data at 757 nm from 2015 to 2017 was used in this study. The footprint size of this dataset was 1.3 km×2.25 km at nadir. Different observation modes of OCO-2 products may cause angular variations of SIF, which impacts relationships between SIF and vegetation production [55]. Following Zhang et al. [70], only OCO-2 SIF data with view zenith angle (VZA) below 20 degrees were used for constraining bidirectional influences while keeping as much available data as possible. The OCO-2 SIF data was classified into soybean and maize groups according to crop type maps and the vertex coordinate of each footprint.
The 16-day composite MODIS EVI product (MOD13A1, collection 6) at a spatial resolution of 500 m from 2015 to 2017 was also used here for assessing the ability OCO-2 SIF to estimate crop yield. We also used a yearly net primary productivity (NPP) product [71] with a spatial resolution of 30 m for the transfer learning as described in Section 2.3.2. The NPP data were produced for the conterminous United States (CONUS) using the MODIS MOD17 algorithm with Landsat surface reflectance, meteorological, and land cover data.

Climate Data
Climate data from 2015 to 2017 were retrieved from the DAYMET database, which provides gridded weather parameters for North America at a spatial resolution of 1 km [68]. The data used include daily precipitation (P), maximum air temperature (T), and water vapor pressure. Vapor pressure deficit (VPD) was further calculated using VP and T [72]. P, T, and VPD were selected according to the previous study in a similar study region [46]. The monthly means of T, P, and VPD were calculated from corresponding daily values for further analysis.

Crop Classification and Yield Data
We used the United States Department of Agriculture (USDA) National Agricultural Statistics Service (NASS) Crop Data Layer (CDL) data to identify areas of corn and soybeans from 2015 to 2017. The CDL data are annual land cover classification products with a spatial resolution of 30 m. The data of county-level soybean and maize yields were retrieved from the USDA Quick Statistic Database.
The CDL data in each year from 2015 to 2017 were employed to determine whether a MODIS EVI pixel and an OCO-2 footprint belong to the soybean or maize category. The area ratios of soybean and maize within all pixels and footprints were calculated. Only pixels and footprints with the area ratio of soybean or maize above 60% were selected for further analysis. A threshold of 60% was used for balancing the tradeoff between the relative purity of the pixels and the availability of the data. A similar value was employed in a previous study to identify vegetation types covered by OCO-2 footprints [62]. For a county, the three-year mean EVI of soybean and maize was calculated for each month in the growing season (from May to September) according to values of selected individual pixels, i.e., where is the three-mean of EVI in month m; EVIy,m,p is the EVI value of the p th pixel in month m of the y th year; Ny is the number of pixels belong to corn or soybean in the y th year.
OCO-2 SIF data was spatially discrete and there are gaps between nearby swaths (Figure 1). Within a county, SIF might be observed on different adjacent days for a given period. The variation of SIF caused by the difference in observing dates must be taken into account. Therefore, the countylevel three-year SIF means of soybean and maize were generated in the following way. All SIF observations of soybean and maize in all three years within each county were lumped together to fit phenology curves of individual crop types (section 2.3.1) since observations of SIF in one year were very limited for a county. Two respective phenology curves of soybean and maize were constructed for each county. With the fitted curves, daily SIF during the growing season was simulated for each county. With simulated daily SIF, monthly means and maxima of SIF were calculated.

SIF Phenology Fitting
A double-sigmoidal function was used to fit the seasonal cycle SIF representing photosynthetic status signals for each county [73]. It has the advantage of retrieving vegetation phenology with noisy data [74,75], such as OCO-2 SIF. This equation is as follows: where y represents the SIF value, x is a given day of the year (DOY); ɑ1 represents the SIF value in the winter; ɑ2 and ɑ3 represent values at the spring and autumn plateau, respectively; b1 and b2 are the DOY mid-points of transitions for spring greenup and autumn browndown, respectively; and d1 and d2 are the corresponding slope coefficients of these transitions. Parameters (ɑ1-ɑ3, b1, b2, d1, d2) on the right of Equation (2) were fitted using three-year mean SIF in different months by a genetic algorithm [73]. Only counties with seasonal curves of SIF well fitted (R 2 > 0.8) were used for further analysis. Finally, there were a total of 218 and 201 counties with SIF data available for maize and soybean, respectively. For all of these available counties, three-year means of EVI and climate data were accordingly calculated for all individual months from May to September.

Model Construction and Validation
The machine learning method adopted here is DNN, one of the typical deep learning algorithms. The DNN has the advantage of learning complex, non-linear relationships and extracting intricate information effectively through the process of transforming raw inputs (satellite and climate data here) gradually towards a higher level of information [34,76,77], e.g., crop yields in this case. DNN used in this study is a feed-forward neural network (NN) with architectures of multiple fully connected layers. Layers in NN can be categorized into the input, hidden, and output layers, and each layer is composed of a certain number of neurons. The number of input neurons is equal to the number of inputs (three-year means of monthly EVI, SIF, and climate variables).
DNN has a large number of hyper-parameters and weights, which require sufficient training data to prevent overfitting and to achieve good model performance. However, due to the low data availability of OCO-2 SIF, it is difficult to collect adequate training, validation, and testing to optimize hyperparameters while testing the trained DNN unbiasedly. Transfer learning (TL) between similar tasks, therefore, becomes desirable for overcoming this problem. Theoretically, this technique enables the leveraging of data and knowledge from a source domain to facilitate modeling in a related target domain, and thus alleviates difficulties and efforts of data collecting and model building [78]. Additionally, the TL technique, unlike the traditional machine learning methods, does not require the training and targeted data in the same feature space and have the same distribution [79]. Based on the rationale of TL, we selected NPP, EVI, and climate data as the source dataset considering the close link between crop NPP and yield. We then built DNNs by transferring the knowledge of employing EVI and climate data to predict NPP to our targeted task of using a similar dataset to predict crop yields. A detailed description of the transferring process is as follows: before the construction and validation of DNN models, all variables were standardized into features with zero mean and unit variance, including yields and NPP of maize and soybean (dependent variables), and three-year means of monthly SIF, EVI, and climate data from May to September (independent variables). Then, we pre-trained the structure of DNN using county-level mean annual NPP and monthly mean EVI from May to September in years 2013 and 2014 to determine hyper-parameters, such as the number of hidden layers and neurons in each layer (Appendix A). The trained network structure is shown in Figure S3. With the DNN architecture obtained above, the weights in each layer were optimized with the training subsets of the targeted dataset (the three-year-mean crop yield, SIF, EVI, and climate data) for corn and soybean, respectively.
The repeated five-fold cross-validation (CV) technique [80] was used to separate the targeted dataset into training and testing groups. In each iteration, the targeted dataset was first shuffled randomly and then split into five equal-sized subsets, of which four subsets of data were used for training the model and the remaining one subset was used for testing the model. The training and testing processes were repeated until all five subsets of data had been used for model tests. All estimated corn and soybean yields in all five test processes were combined to produce a complete validation dataset for this iteration. This whole process of five-fold CV was then iterated 100 times. The validation scores of the 100 runs were averaged for further analysis.
The performance of models was evaluated with statistical metrics, including R 2 , mean absolute error (MAE), mean absolute percentage error (MAPE), and percentage error (PE): where is the crop yield estimated by the model, is the three-year-mean crop yield calculated based on measured data from USDA, indicates a given county, and n is the number of available records.

Importance Determination of Different Variables
The feature importance (FI) analysis was employed to identify the importance of different variables in crop yield prediction. The 'feature' here refers to model predictors, including SIF, EVI, and climate variables. The FI technique randomly permutates the values of one feature while keeping the values of remaining features unchanged and then evaluates the relative importance of the aimed feature by measuring the reduction of model performance [81]. The rationale behind this algorithm is that by shuffling the feature values randomly, meaningful information of the feature is replaced with random noise. As a result, the association between feature and response variables is broken, leading to the decrease of the variance explained by the model (R 2 ) [82,83]. In the current study, we calculate FI scores by averaging the results from the 100 repeated experiments, and the FI score for a feature is defined as follows: where is the R 2 of the original model without permutation and is the value of the model after permutating the feature variable; M is the number of times the feature is shuffled, which aims to increase the robustness of the method [81].

Relationships Between Crop Yields and Different Variables
The relationships of crop yields with monthly SIF, EVI, and climate variables were quantified using the Pearson product-moment correlation coefficients (r) (Figure 3). Monthly VPD generally had negative correlations with both corn and soybean yields, especially in the peak growing season (June, July, and August) (p<0.05) (Table S5). Precipitation in all months was positively correlated with both corn and soybean yields. Correlations of monthly means and maxima of SIF and EVI with corn and soybean yields sharply increased from May to June and then decreased from August to October. In September and October, SIF showed stronger positive correlations with the yields than EVI. In most of the months, both monthly maxima of SIF and EVI performed better in indicating the yields than the corresponding monthly means. Therefore, monthly maxima of SIF and EVI were used in further analysis.

Performances of DNNs for Corn and Soybean Yield Prediction
County-level yields of corn and soybean, monthly means of climate (T, P, and VPD) and maxima of EVI and SIF from May to September were used to train DNN models. Table 1 shows the performances of DNN models with different combinations of inputs. Models with inputs of SIF, EVI, SIF plus EVI, and Climate have moderate to high coefficients of determination, indicating the usefulness of satellite and climate information in predicting corn and soybean yields. Among models only with EVI, SIF, and climate as inputs, the models with climate predictors performed the best, explaining 76% and 82% variations of county-level corn and soybean yields, respectively. Models only with monthly maximum SIF as inputs performed the poorest. The integration of SIF and EVI could improve the performance of yield estimating models that only used either SIF or EVI data as predictors. In addition, the combination of satellite and climate data significantly improved model performances. Specifically, models with both EVI and climate as inputs outperformed the climatebased models, increasing R 2 from 0.76 to 0.86 for corn and from 0.82 to 0.87 for soybean. The combination of climate information with SIF also improved the estimates of the yields, with R 2 increased from 0.76 to 0.80 for corn and from 0.82 to 0.83 for soybean. These results suggest that OCO-2 SIF could make some contributions to the yield estimate, but less than MODIS EVI. The models with all SIF, EVI, and climate variables as inputs performed similar to the models with EVI and climate as inputs, implying that OCO-2 SIF provided less useful information for capturing the spatial variations of corn and soybean yields in the study area than EVI. The analysis of variance (ANOVA) and the Tukey Honest Significant Difference (HSD) tests [84] were conducted to investigate whether the differences among model performances in the 100 repeated experiments are significant (Table S2 and Table S3). The ANOVA tests for both corn and soybean showed significant differences among performances of different models (p < 0.001). Additionally, the Tukey HSD tests demonstrated that most of the performance differences between paired models were significant, except for the differences between 'EVI+Climate' and 'SIF+EVI+Climate' models for both corn and soybean, and between 'Climate' and 'SIF+EVI' for corn. These results further confirmed that the addition of SIF data into the 'EVI+Climate' models did not improve the estimates of both corn and soybean yields. However, there are significant differences between the 'EVI' and the 'SIF+EVI' models, and between the 'Climate' and the 'SIF+Climate' models. This indicated that the inclusion of SIF could improve the 'Climate' model and 'EVI' model for estimating yields of corn and soybean.
The performances and stability of DNNs were compared by summarizing the repeated CV results. As shown in Figure 4, the 'Climate' model had lower minimum and first quantile values of R 2 than those of the 'SIF+Climate' and 'EVI+Climate' models for both corn and soybean. This means that, even for the worst cases, the models with both satellite and climate predictors outperformed the models with only climate data as inputs. A similar finding was also observed between the 'SIF+EVI' and 'EVI' models for soybean, and the former performed better than the latter under the worst and other conditions. The interquartile range (IQR) and overall range of R 2 of 'EVI+Climate' and 'SIF+EVI+Climate' models were noticeably smaller than those of the 'Climate' models for both corn and soybean (Table S4), indicating that the overall performance of 'EVI+Climate' and 'SIF+EVI+Climate' models varied at a lower degree among different experiments than the 'Climate' model. This result demonstrated that the addition of EVI and SIF data into the climate-based models enhanced the stability of DNN performances and the reliability of yield estimates under different conditions.  Figure 5 and Figure 6 show the spatial distribution of PE for different models. The spatial patterns of PE for both corn and soybean were similar among the 'EVI+Climate,' 'SIF+Climate,' and 'SIF+EVI+Climate' models. In particular, those models performed similarly well in the central part of the Corn Belt (e.g., Iowa). The magnitudes of PE were mostly smaller than 2%. In contrast, PE was higher in the northern part of the Midwest states (i.e., North Dakota, northern South Dakota, and northwestern Minnesota), which indicates that models tended to underestimate yields of soybean in these regions (Figure 5b, d and f). The underestimation of corn and soybean yields by 'EVI', 'SIF', and 'SIF + EVI' models were also significant in the northwestern part of the study area (Figure 6b, d, and f). In contrast, the PE values of models with only climate variables as inputs showed less spatial variability (Figure 6h).

Spatial Differences in Performances of the DNN Models
For corn, the PE values of models based only on SIF, EVI, and climate inputs ranged from -20.25% to 47.03%, from -26.11% to 46.09%, and from 15.68% to 30.85% (Table 2). It was also found that the combination of remote sensing and climate data made the range of PE much smaller, i.e., -12.46% to 21.90 for the model based on SIF and climate, -12.28% to 15.93% for the model based on EVI and climate, and -11.44% to 17.83 for the model with SIF, EVI, and climate data as inputs. As for soybean, the PE ranges by individual models are similar to those of corn. In addition, the DNN models performed well at the county level for both crop types. For above 70% of counties, PE values of 'EVI+Climate' and 'SIF+EVI+Climate' models ranged from -3.48% to 3.75% and from -3.55% to 3.39% for corn. The corresponding values were from -4.37% to 4.65% and from -3.63% to 4.27% for soybean ( Table 2).  and right are for the soybean. The figure includes results from models that input both SIF and climate data (a and b), input both EVI and climate data (c and d), and input SIF, EVI, and climate data (e and f). In each plot, the density plot in the upright corner presents the overall distribution of PE. The red dashed lines and 'M' in the density plots represent the mean values of PE, and the black ones are the referenced lines located at zero. Positive values mean that estimated yields were lower than census data, vice versa.

Importance of Different Variables
The feature importance (FI) scores of predictors in the 'EVI + SIF + Climate' models were calculated for examining the importance of individual predictors in estimating corn and soybean yields. The 'EVI + SIF + Climate' models were used because they included the most comprehensive inputs relative to other models using only satellite or climate predictors. Figure 7 shows the importance of the variables in predicting corn and soybean yields, with a higher FI value indicating a greater contribution to the model performances. It should also be noted that the FI analysis is modelspecific because different variables might have dissimilar functions in different DNNs, and the FI results for other models were presented in Figure S4-6.
In the 'EVI+SIF+Climate' model, both satellite and climate predictors in the peak growing months (from June to August), which is corresponding with the grain filling and reproductive periods in the Midwestern U.S. [85], were essential in estimating crop yields. Specifically, EVI in July was the most important variable determining yields of both corn and soybean. In this month, corn and soybean here started the reproductive stage [45,86], and EVI effectively captures the growing status of the crop. EVI in August is the second important factor for predicting soybean yield and the third important remote sensing factor for predicting corn yield in the Midwest states, which was also confirmed by a previous study [46]. They declared that the inclusion of EVI in August noticeably reduced the uncertainty of yield prediction. As for climate variables, VPD in July and temperature in June acted as the important determinants for yields of corn and soybean, respectively. VPD has a strong impact on the stomatal conductance of crops, which affects carbon assimilation [7,87,88]. As shown in Figure 3, VPD in July had a significant negative correlation with the corn yield. Temperature determines VPD and influences photosynthetic activities, and thus constrains biomass production and influences crop yield [89][90][91]. In the study region, the temperature in June played an important role in determining the yield of soybean, just following EVI in July and August (Figures 3  and 7). SIF was generally less important for estimating corn and soybean yields in comparison with EVI and climate variables. The information carried by SIF into the model might be overshadowed by EVI. As shown in Table 1, 'EVI+SIF+Climate' and 'EVI+Climate' models performed very similarly.

4.Discussion
SIF is tightly linked with vegetation photosynthesis at various temporal scales. Theoretically, SIF has the potential to be an effective predictor of crop yields, because it is associated with electron transport in light reaction of photosynthesis, and therefore, closely related to GPP of plants [48,49,64]. However, the spatial resolution of most available satellite SIF is relatively low (40 km×80 km for GOME-2 SIF and 30 km×60 for SCIAMACHY SIF), and it is difficult to use such low-resolution satellite SIF to estimate crop yields at small scales, such as the county-level. The spatial resolution of OCO-2 SIF is much higher than those of GOME-2 and SCIAMACHY. However, the swath width of OCO-2 is only 10.3 km. There are gaps around 100 km between two adjacent swaths [92]. The spatial and temporal discontinuities of OCO-2 SIF hinder the application of these data in estimating crop yields at regional scales.
In this study, three years of OCO-2 SIF data in different months were used to fit the phonology curves of corn and soybean for individual counties. With the fitted curves, monthly SIF means and maxima were estimated for corn and soybean in different counties. These monthly SIF values in the peak growing season were significantly positively correlated with three-year mean county-level corn and soybean yields (Figure 3), confirming that SIF is an effective indicator of crop yields owing to its direct linkage with photosynthesis. However, compared to EVI, our processed SIF had relatively weaker correlations with corn and soybean yields ( Figure 3) and poor ability to capture the spatial variations of yields (Table 1). This is due to the fact that the county-level EVI used here was calculated from spatially and temporally continuous direct observations by the satellite, while the SIF values used were estimated from the phenology curves which were fitted from limited direct observations in three years. The numbers of SIF observations were uneven for a county, and locations of footprints might differ in different years. The discrepancy between estimated and real county-level mean SIF would definitely weaken the effectiveness of OCO-2 SIF data in predicting crop yields. Therefore, the development of spatially and temporally continuous high-resolution SIF datasets is required for improving estimates of crop yields.
In the study area, climate-based DNN models overall performed better than remote sensing-based DNN models in capturing spatial variations of county-level corn and soybean yields, possibly due to the strong controls of VPD and precipitation on yields (Figure 3). In addition, the departures of yields estimated by climate-based DNN models from census data showed less spatial variations relative to those form SIF and EVI-based DNN models ( Figure 6). These results indicated that climate information played a more important role in controlling the spatial variations of corn and soybean yields in this region. EVI was an effective indicator of crop biomass, which is tightly linked with final yields. Therefore, the performance of EVI-based DNN models was close to that of climate-based DNN models (Table 2 and Figure 4). However, similar to SIF-based models, EVI-based models performed poorer than climate-based models in counties with low area fractions of corn and soybean ( Figure 6 and Figure S2), indicating that there were impacts of land cover heterogeneity on crop estimates by remote sensing data.
The combination of different sources of data might improve crop yield estimation (Table 1). Although, when compared with other model variables, our estimated county-level mean SIF showed relatively poor ability to capture spatial variations of corn and soybean yields, it could act as effective supplement information of EVI and climate for better estimating crop yields. DNN models inputting both SIF and EVI data performed very close to climate-based DNN models. Additionally, DNN models integrating EVI and climate data outperformed either EVI or climate-based models (Figures 4, 5, and 6; Tables 1 and 2). However, the addition of data from the same source could not constrain yield estimate. From the statistical point of view, DNN models with SIF, EVI, and climate data as inputs performed even slightly poorer than models with EVI and climate ( Table 1). The spatial patterns of discrepancies between yields estimated by two models and census data were also similar ( Figure 5).
Nevertheless, satellite SIF data is useful for estimating crop yields even if current SIF products have some limitations. The launching of satellite with finer resolution and broader spatial coverage than OCO-2, such as FLEX [93] and GeoCARB [94], would pave the way for further testing to what degree SIF might contribute to crop yields estimation from space. Currently, one possible solution for estimating crop yield using SIF data is to downscale satellite SIF data. For example, several studies have proposed methods to produce spatially continuous high spatial resolution SIF datasets using OCO-2 SIF measurements [95][96][97]. It is worthy of investigating the applicability of such downscaled SIF in estimating crop yields.
Furthermore, it should be noted that when predicting crop yield empirically, especially when employing sophisticated machine learning techniques such as DNN, it would be beneficial to collect more training data. In addition, information other than EVI or SIF could also be useful in predicting corn and soybean yield and boosting the model performance as additional predictors. For example, vegetation optical depth (VOD) derived from the satellite passive microwave bands could indicate canopy biomass and water content condition of the plants [35,98,99]; evaporative stress index (ESI) retrieved from satellite thermal bands and leaf area index is desirable for assessing the crop moisture status as an indicator of agricultural drought [100,101]. Furthermore, it is expected that with more data and more comprehensive predictors, the model could extract information more effectively and better discover non-linear relationships between dependent and independent variables. Therefore, our future work would further test the usefulness of SIF signals in crop yield predicting when more data are available. We would also employ multiple sources of data for crop yield predicting via the DNN technique, aiming to build a high-performance crop yield forecasting system at the regional level.

5.Conclusion
In this study, we proposed a method to use the OCO-2 SIF product to predict yields of corn and soybean in the U.S. Corn Belt. The ability of OCO-2 SIF products to estimate county-level crop yields in the study area was compared with those of MODIS EVI and conventional climate data. The following conclusions could be drawn from this study.
(1) Monthly mean and maximum SIF during the peak growing season showed similar significant correlations with corn and soybean yields to MODIS EVI. The DNNs with SIF as predictors were able to capture spatial variations of corn and soybean yields but performed much poorer than those based on MODIS EVI and climate variables. The performance of SIF and MODIS EVI-based DNNs varied with the areal dominance of crops. In areas with low areal fraction of crops, the SIF and EVI-based DNNs tended to underestimate corn and soybean yields noticeably. The performance of DNNs with climate variables as predictors exhibited less spatial variability.
(2) The combination of remote sensing and climate data improved the estimation of crop yields. Our processed county-level SIF data could provide useful supplementary information to MODIS EVI and climatic variables in estimating crop yields. DNNs of EVI and climatic variables in conjunction with SIF obliviously outperformed DNNs with only EVI or climatic variables as predictors. The integration of MODIS EVI and climatic variables was able to improve the estimating accuracy of crop yields significantly. Further addition of SIF into models with MODIS EVI and climatic variables only marginally enhanced crop yields prediction.
(3) Feature importance analysis indicated that MODIS EVI and climate predictors (e.g., VPD and temperature) during the peak growing season (from June to August) played critical roles in predicting yields of corn and soybean in the Midwestern 12 states in the U.S. for the DNNs with all SIF, EVI, and climatic variables as predictors. By contrast, the importance of SIF was much smaller than MODIS EVI and climatic variables in those DNNs.
Supplementary Materials: The following are available online at www.mdpi.com/2072-4292/12/7/1111/s1, Figure S1: Corn and soybean harvested area fraction and production fraction of the 12 Midwest states in the U.S. from 2015 to 2017. The data were compiled from USDA NASS crop data layer (CDL) and crop yield datasets. Figure S2: Corn and soybean growing area fraction in each county. The fractions were three-year-mean values (from 2015 to 2017) calculated using USDA NASS CDL. APE (absolute percentage error) in the barplots were calculated by averaging county-level absolute PE values of 'SIF', 'EVI', and 'Climate' models. In the groups with area fraction equal to 5% -15%, APE values of both 'SIF' and 'EVI' models were significantly lower than the 'Climate' models for corn and sobyean. Figure S3: DNN structure employed in this study. The neural network consisted of one input layer, three hidden layers, and one output layer. (n = the number of input variables; m = the number of records). Figure S4: Feature importance (FI) scores of predictors of the 'SIF+EVI' model for corn (left panels) and soybean (right panels). Values of the x-axis are feature importance scores, and the variable with larger FI score carries more importance in predicting the crop yield. Figure S5: Feature importance (FI) scores of predictors of the 'EVI + Climate' model for corn (left panels) and soybean (right panels). Values of the x-axis are feature importance scores, and the variable with larger FI score carries more importance in predicting the crop yield. (P = precipitation; T = temperature; VPD = vapor pressure deficit). Figure S6: Feature importance (FI) scores of predictors of the 'SIF + Climate' model for corn (left panels) and soybean (right panels). Values of the x-axis are feature importance scores, and the variable with larger FI score carries more importance in predicting the crop yield. (P = precipitation; T = temperature; VPD = vapor pressure deficit). Table S1: Trained and tested R 2 of different model architectures. The best DNN has a hidden layer of three with 35 neurons in each layer. Table S1: Summary of Tukey HSD test between adjusted R 2 values of different models in 100 experiments for corn. 'Group1' and 'Group2' are the targeting pair of models for comparison. 'MeanDiff' is the difference between group means. 'Lower' and 'Upper' are the lower and upper boundaries of confidence intervals for the pairwise mean differences. 'Reject' is TRUE means that there is a significant difference (at 95% confidence level) between two models; FALSE means that the null hypothesis is true that the averaged model performances measured by R 2 are equal. Table S3: Summary of Tukey HSD test between adjusted R 2 values of different models in 100 experiments for soybean. 'Group1' and 'Group2' are the targeting pair of models for comparison. 'MeanDiff' is the difference between group means. 'Lower' and 'Upper' are the lower and upper confidence interval boundaries for the pairwise mean differences. 'Reject' is TRUE means that there is a significant difference (at 95% confidence level) between two models; FALSE means that the null hypothesis is true that the averaged model performances measured by R 2 are equal. Table S4: Summary of R 2 values of 100 experiments of different models. Q1, Q2, and Q3 represent the first, second, and third quantiles separately; IQR is the interquartile range; max = Q1-1.5*IQR, min = Q3+1.5*IQR; range = max-min. Table S5: P values of correlation relationships between variables and crop yield. The shaded area represent the predictors that were included in DNN models. Bold numbers are the model predictors of which the p values are less than 0.05.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A: Pre-train DNN and the Auxiliary Dataset
Considering that transfer learning enables leveraging knowledge from similar tasks, we designed an auxiliary dataset in order to pre-train the DNN to obtain the optimal neural network architecture that is suitable for our targeted task, i.e. predicting crop yield. We used the net primary production (NPP) product derived from the Landsat [71] as an approximate of the crop yield, the dependent variable in our targeted task; and the independent variables include EVI, VPD, precipitation, and temperature, which are the same as our targeted dataset. The auxiliary dataset from 2013 to 2015 was processed in the GEE using the same method described in section 2.2.3 to obtain county-level mean values. Then we split the auxiliary dataset into two parts. Specifically, the first part of the data was from 2013 to 2014, and the second set, or the testing set, included data in 2015. Next, the 5-fold CV was employed on the first set to select training and validation data and to optimize hyperparameters. The hyperparameters obtained mainly included the number of layers of the neuron network, the number of neurons in each layer, the activation function, the dropout rate, and the optimize algorithm. The best hyperparameters were selected based on the validation dataset by choosing the ones with the highest averaged R2. Last, the DNN with best hyperparameters was tested on the second set. Part of the testing results for model architecture selection was presented in Table S1, and the model structure was shown in Figure S3. It should also be noted that the weights from the pre-trained model of which the NN architecture was selected are not used for initialization. Instead, the DNN models in Section 3.2 were all randomly initialized in order to prevent the information from EVI as the pre-trained dataset contaminate the 'SIF' models and to provide a relatively fair comparisons among models that input EVI and the ones that input SIF data.