A Satellite-Based Method for National Winter Wheat Yield Estimating in China

Abstract

Satellite-based models have tremendous potential for monitoring crop production because satellite data can provide temporally and spatially continuous crop growth information at large scale. This study used a satellite-based vegetation production model (i.e., eddy covariance light use efficiency, EC-LUE) to estimate national winter wheat gross primary production, and then combined this model with the harvest index (ratio of aboveground biomass to yield) to convert the estimated winter wheat production to yield. Specifically, considering the spatial differences of the harvest index, we used a cross-validation method to invert the harvest index of winter wheat among counties, municipalities and provinces. Using the field-surveyed and statistical yield data, we evaluated the model performance, and found the model could explain more than 50% of the spatial variations of the yield both in field-surveyed regions and most administrative units. Overall, the mean absolute percentage errors of the yield are less than 20% in most counties, municipalities and provinces, and the mean absolute percentage errors for the production of winter wheat at the national scale is 4.06%. This study demonstrates that a satellite-based model is an alternative method for crop yield estimation on a larger scale.

1. Introduction

As one of the largest producers of wheat, China plays a dominant role in retaining the global balance of demand and supply of wheat [1,2]. In 2018, China produced 131.44 million tons of wheat, which accounted for around 17.9% of the global wheat production [3,4]. Winter wheat makes up approximately 95.6% of the wheat production in China according to the 2019 statistical data (https://data.stats.gov.cn, accessed on 16 October 2021). The growing season of winter wheat coincides with the dry season (November to June), suggesting that there is a high possibility that wheat experiences dry events [5]. Therefore, reliable estimations of winter wheat production in China on a national scale are extremely important for governments when making decisions regarding food security.

A variety of methods have been used for estimating crop yield: statistical models, machine learning algorithms, crop growth modeling and light use efficiency (LUE)-based models. Statistical models rely on the regressions between crop yield and various variables, such as the normalized difference vegetation index (NDVI) [6,7], the enhanced vegetation index (EVI) [8], and land surface temperature [9,10]. Statistical models usually show good performance for crop yield estimation in a given study area, but can hardly be applied to large areas owing to the limited spatial generalization [11,12]. Machine learning algorithms do not require the presupposed relationships between yield and predictor variables [13,14], and usually show good performance [15]. For example, a recent study highlighted that machine learning algorithms could predict the yield at the county level about four months prior to harvest [16]. However, the performance of machine learning algorithms usually needs large volume of observations for model training, and the models are specific for a given region and period and unavailable for direct support to other areas and phases [17,18]. Crop growth modeling was developed to represent the critical processes of crop growth as well as the responses of crops to the environment [19,20,21]. Due to the complexity of crop growth and existing knowledge gaps, there are still large challenges for large-scale applications of crop growth modeling [22,23,24]. For example, crop growth modeling requires a large set of model parameters to simulate the physiological and growth processes of crops, which makes them difficult to implement for large-scale applications [25,26].

Satellite-based methods have been widely used to simulate crop production over large areas, benefitting from temporally and spatially continuous crop growth information derived from satellite data [27,28]. LUE models, which are based on satellite data, are a powerful tool for quantifying crop yield on a large scale [29,30]. LUE models are designed to simulate vegetation gross primary production (GPP) based on the assumption that GPP is directly dependent on the absorbed photosynthetically active radiation (APAR) through LUE [31,32]. Previous studies have validated the performance of the LUE model in simulating GPP of major ecosystem types at regional and global scales [33,34,35,36,37]. Then, LUE models were further developed to simulate crop yield by calculating the ratio of GPP transferred to harvested crop organs [38]. Several studies have employed LUE models to simulate the yield of various crops and their spatial and temporal variations at the regional scale [27,28,39]. For example, He et al. [40] applied a LUE model to produce 30 m spatial resolution GPP of seven crop types in Montana and estimate their yield, and indicated that the estimated production was consistent with reported production at the county scale. A recent study also showed that a satellite-based LUE model can capably simulate the temporal and spatial patterns of winter wheat yield in Kansas [41], which can be used for future crop yield estimation on a larger scale.

Notably, it is still a huge challenge to estimate crop yield over large regions or at a national scale [42,43], because the crop yield is highly dependent on crop variety and management practices [44,45], which show significant heterogeneity [46]. López-Lozano et al. [47] found strong correlations between crop yield and accumulated fAPAR in Europe, but these correlations largely differed with crop types. Based on MODIS-GPP production (MOD17), Reeves et al. [48] revealed that the estimated GPP and reported wheat yield had a weak correlation at county scale and climate district in the United States and emphasized the importance of using spatial-specific harvest index. Therefore, it is still a challenge to estimate crop yield with large spatial heterogeneity, especially on a national scale. Here, this study aims to examine the ability of a satellite-based LUE model (i.e., EC-LUE) [49,50] to quantify the yield and production of winter wheat in China. The overarching goals of this study are to (1) examine the model performance for simulating crop GPP based on eddy covariance (EC) measurements, (2) assess the accuracy of EC-LUE for simulating winter wheat yield at field, county, municipal and province scales, and (3) investigate the ability to quantify the production of winter wheat at the national scale.

2. Materials and Methods

2.1. Study Area

The study area covers 11 provinces: Anhui (AH), Gansu (GS), Hebei (HB), Henan (HN), Hubei (HuB), Jiangsu (JS), Shandong (SD), Shanxi (SAX), Shanxi (SX), Sichuan (SC) and Xinjiang (XJ). This area includes 146 municipalities and 789 counties with statistical yield information (Figure 1). The winter wheat planting area and production in our study region account for 96% and 98% of the total planting area and production in China, respectively (China Rural Statistical Yearbook 2020). Generally, winter wheat is sown in early-mid October and harvested in early-mid June of the following year.

2.2. Data

2.2.1. Model Forcing Data

In this study, we derived the NDVI for 2016–2020 from the Enhanced Thematic Mapper Plus (ETM+) sensor onboard Landsat 7, the Operational Land Imager (OLI) sensor onboard Landsat 8 (United States Geological Survey, https://www.usgs.gov/, accessed on 16 October 2021) and Multispectral Instrument (MSI) sensor onboard Sentinel 2 (European Space Agency, https://github.com/senbox-org, accessed on 16 October 2021), which were processed for atmospheric correction. In order to increase the observation frequency of satellite as much as possible and reduce the impact of clouds, we firstly resampled the NDVI of Sentinel 2 with a nearest neighbor method to 30 m to keep a consistent spatial resolution with the Landsat image. Then, we obtained the NDVI values of all valid pixels within 16 days after cloud removal, and acquired the maximum value for 16 days synthetized with a spatial resolution of 30 m. The above operations were run on the Google Earth Engine platform. In addition, we used a linear interpolation method to fill the missing values. Furthermore, a Savitzky–Golay filter [51], which has the ability to retain the vegetation signal information, was used to smooth the NDVI series.

The meteorological forcing data of the EC-LUE model, including photosynthetically active radiation (PAR), air temperature (T) and vapor pressure deficit (VPD), were derived from ERA-Interim of the European Centre for Medium-Range Weather Forecasts (ECMWF) (https://cds.climate.copernicus.eu/, accessed on 16 October 2021). VPD was calculated based on T and dew-point temperatures from Yuan et al. [50]. The original spatial resolution of ECWMF was 0.1 × 0.1 degree, and in order to keep the same spatial resolution with NDVI, the PAR, T and VPD were also resampled to 30 m based on the nearest neighbor method and the corresponding 16-day mean values were also obtained.

2.2.2. Eddy Covariance Measurement

Three eddy covariance (EC) measurements, located in Luancheng (LC, 114.6928° E, 37.8894° N) [52], Guantao (GT, 115.1274° E, 36.515° N) [53] and Shouxian (SX, 116.47° E, 32.36° N), were used to examine and calibrate the EC-LUE model (Figure 1). These observed sites are croplands with a wheat-maize (LC and GT) and wheat-rice (SX) rotation. In this study, the nighttime flux partitioning method [54] was applied to calculate the daily average GPP of winter wheat in the whole growing season. Given the uniqueness of the agroecosystem, we used the Markov chain Monte Carlo (MCMC) method to calibrate the value of maximum light use efficiency (ε_max) by comparing estimated GPP based on the EC-LUE model to the calculated GPP from EC measurements. The details of this method can be found in previous publications [55,56].

2.2.3. Winter-Wheat-Related Data

In this study, the 30 m winter wheat planting area maps in China between 2016 to 2018, obtained from Dong et al. [57] (https://doi.org/10.6084/m9.figshare.12003990, accessed on 16 October 2021), were used to estimate the yield of winter wheat. These maps were acquired with a phenology-based planting area identification method, which can realize the early winter wheat identification in April with an overall accuracy of 89.88%.

We collected field survey winter wheat yield data reported by farmers at 30 sites in Henan, Hebei and Shandong provinces from 2019 to 2020 (yellow point in Figure 1), and each site only had one year of data. We used these data to evaluate the performance of the EC-LUE model at the surveyed fields. Specifically, we randomly selected 60% of the surveyed yield data to calibrate the harvest index and then we used the calibrated HI to estimate yield at the remaining sites. Moreover, to further evaluate the model performance in estimating winter wheat yield at province, municipal, and county scales, we collected the statistical winter wheat yield data from 2016 to 2018 from the Statistical Yearbook of each province, and acquired statistical yield data of 11 provinces, 146 municipalities and 789 counties.

2.3. Method of Winter Wheat Yield Estimation

This study used a satellite-based light use efficiency (LUE) model (i.e., eddy covariance light use efficiency, EC-LUE) to estimate the yield and production of winter wheat over the study area. The method first estimates the vegetation gross primary production (GPP) and then converts GPP estimates to crop yield based on the harvest index [38,41,58]:

$$\mathit{Yield}=\mathit{GPP}\times \mathit{AR}\times \frac{1}{1+\mathit{RS}}\times \frac{1}{1-\mathit{MC}}\times \mathit{HI}$$

(1)

where Yield is the estimated winter wheat yield. GPP represents the accumulated GPP of winter wheat throughout the growing season. AR is the remaining component after considering the consumption of autotrophic respiration, which is set to 0.66 for winter wheat [52]. RS represents the ratio of root to aboveground biomass, with a value of 0.2 [58]. MC indicates the moisture content of grain, which is set to 0.11 for winter wheat [59]. HI refers to the harvest index, that is, the proportion of aboveground biomass converted into crop yield [60,61].

The EC-LUE model, developed by Yuan et al. [49,50], was used to simulate the GPP of winter wheat in this study area. An advantage of the EC-LUE model is its ability to simulate GPP over large areas as the potential LUE is independent of various land cover types [62]. The model is driven by the normalized difference vegetation index (NDVI), photosynthetically active radiation (PAR), air temperature (T) and vapor pressure deficit (VPD). The details are as follows:

$${\mathit{GPP}=\mathit{PAR}\times \mathit{fPAR}\times \epsilon }_{\mathit{max}}{\times \min (T}_s{,W}_s)$$

(2)

$$\mathit{fPAR}=1.24\times \mathit{NDVI}-0.168$$

(3)

$$T_S=\frac{\left({T-T}_{\mathit{min}}\right){\times (T-T}_{\mathit{max}})}{\left({T-T}_{\mathit{min}}\right)\times \left({T-T}_{\mathit{max}}\right) {-(T-T}_{\mathit{opt}}{)}^2}$$

(4)

$$W_S=\frac{{\mathit{VPD}}_{\mathit{max}}-\mathit{VPD}}{{\mathit{VPD}}_{\mathit{max}}{-\mathit{VPD}}_{\mathit{min}}}$$

(5)

where fPAR is the part of absorbed PAR. ε_max indicates the potential light use efficiency without environmental stress (2.14 g C m⁻² MJ⁻¹ APAR); min (T_s, W_s) refers to the minimum value between T_s and W_s. The model assumes that the limiting factors (temperature and moisture) follow Liebig’s law, that is, the LUE depends on the most limiting factor. T_min, T_max and T_opt are the minimum, maximum and optimum of air temperature (°C) for vegetation photosynthesis, and their values are 0 °C, 40 °C and 21 °C, respectively. If the air temperature is lower than T_min or higher than T_max, T_s is set to 0 [49]. VPD_min and VPD_max are the minimum and maximum of vapor pressure deficit (Pa), with values of 650 Pa and 4300 Pa, respectively. When VPD > VPD_max, W_s is 0, and when VPD < VPD_min, it is 1 [63].

In this study, we first estimated GPP at the 30-m spatial resolution based on the newest planting area map of winter wheat by Dong et al. [57] (Section 2.2.3). Considering the temporal and spatial heterogeneity, we calculated the HI of each province, municipality, and county based on statistical yield data and estimated GPP. We used statistical data over a 3-year period (2016–2018) and employed the cross-validation method to invert the HI and validate model performance. Specifically, we used the two years of statistical data and estimated GPP to invert the HI in each province, municipality, and county, and used the inverted HI to estimate the yield based on Equation (1) for the remaining year to examine model performance. This procedure was repeated three times and each year was used as an independent validation.

In order to estimate the production of winter wheat at the national scale, we first estimated the winter wheat yield in each province base on the accumulated GPP in the growing season and the inverted harvest index at a province level. Second, we computed the planting area of winter wheat in each province according to the winter wheat map. Finally, we multiplied the estimated yield by planting area to obtain the total production of all provinces in the corresponding year, according to:

$$\mathit{Prod}={{\displaystyle \sum }}_{i=1}^n{(\mathit{Yield}}_i{\times \mathit{Area}}_i)$$

(6)

where Prod represents the total production of all provinces, i refers to the number of provinces (i = 1, 2, 3…11), Yield is the estimated yield and Area is the planting area of winter wheat in the corresponding province.

2.4. Model Evaluation

We evaluated the performance of our method for reproducing GPP, yield and production based on four metrics: the coefficient of determination (R²), root mean square error (RMSE), mean absolute percentage error (MAPE) and the refined Willmott’s index of agreement (d_r).

$$\mathit{RMSE}=\sqrt{ \frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{\left(P_i-O_i\right)}^2}$$

(7)

$$\mathrm{MAPE}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n|\frac{O_i-P_i}{O_i}|\times 100$$

(8)

$$d_r=\{\begin{array}{c}1-\frac{{{\displaystyle \sum }}_{i=1}^n\left|P_i-O_i\right|}{c{{\displaystyle \sum }}_{i=1}^n\left|O_i-\overline{O}\right|} , when{\displaystyle \sum }_{i=1}^n\left|P_i-O_i\right| \le c{\displaystyle \sum }_{i=1}^n\left|O_i-\overline{O}\right|\\ \frac{c{{\displaystyle \sum }}_{i=1}^n\left|O_i-\overline{O}\right|}{{{\displaystyle \sum }}_{i=1}^n\left|P_i-O_i\right|}-1, when{\displaystyle \sum }_{i=1}^n\left|P_i-O_i\right|>c{\displaystyle \sum }_{i=1}^n\left|O_i-\overline{O}\right|\end{array}$$

(9)

where n is the number of administrative areas, P_i and O_i are the estimated and statistical yield at the ith administrative area, respectively. $\overline{O}$ is the mean of the statistical yield over all administrative areas, and c equals 2, d_r ranges from −1 to 1, the greater value indicating better model performance [64].

3. Results

3.1. GPP Simulation

We calibrated the model parameter (i.e., maximum light use efficiency, ε_max) at the LC, GT and SX sites by the MCMC method, with the values equal to 3.43 g C MJ⁻¹, 3.36 g C MJ⁻¹ and 3.38 g C MJ⁻¹, respectively. Specifically, we used the ε_max of SX to validate the model performance for simulating GPP at other two sites. Overall, the simulated GPP from the EC-LUE model agrees well with the calculated GPP based on EC measurements (Figure 2), suggesting that the model is capable of capturing the temporal variations and magnitude of tower-based GPP. For example, at the LC site, which has the longest observations from 2008 to 2016, the coefficient of determination (R²) is 0.84 (Figure 2d). In addition, the d_r are 0.78, 0.8 and 0.85 at the LC, GT and SX sites, respectively.

3.2. Yield and Production Estimation

We first evaluated the model performance in estimating the yield at the site scale, based on field-surveyed yield data at the 30 sites of the Henan, Hebei and Shandong provinces. Specifically, we first calibrated the harvest index (HI, 0.48) and then used the calibrated HI to estimate yield and examined the model performance (Section 2.2.3). The results show that the model is able to explain 51% of the spatial variations in investigated yield (Figure 3). The RMSE, MAPE and d_r are 28.93 g C m⁻² yr⁻¹, 5.85% and 0.64, respectively.

In addition, we examined the model performance for estimating winter wheat yield at municipal and county scales based on the statistical dataset. At the municipal scale, the model can simulate the spatial variations of yield over all 11 provinces (Figure 4). The R² between statistical and estimated yield ranges from 0.47 at HB and XJ to 0.8 at SAX (Figure 4 and Figure 5b), and all slopes of the linear regression are larger than 0.81 (Figure 5a), indicating no obvious systematic estimates errors at the municipal scale. The RMSE ranges from 19.9 g C m⁻² yr⁻¹ to 45.06 g C m⁻² yr⁻¹, the MAPE values are less than 20% and d_r values are greater than 0.5 in most of the investigated provinces (Figure 5c–e), implying a smaller deviation between the statistical and estimated yield. However, there still exists some differences between statistical and estimated yield in several provinces. For example, poor performance is observed in HB, JS and XJ, although the MAPE values are smaller than 20% (Figure 5d), and the R² and d_r are both less than 0.5 (Figure 5b,e).

Overall, at the county scale, the model estimations can simulate the spatial variations of the statistical yield in most provinces (Figure 6), with the R² between statistical and estimated yield varying from 0.37 to 0.82 (Figure 5b and Figure 6), and the RMSE ranging from 23.73 g C m⁻² yr⁻¹ to 46.65 g C m⁻² yr⁻¹ (Figure 5c) (except for the AH, GS and XJ provinces). Moreover, the MAPE values are lower than 20% and d_r values are higher than 0.5 in most provinces (Figure 5d,e), which also shows a lower deviation between statistical and estimated county yield. Although the model produces praiseworthy performance at the county scale, there are large uncertainties over several provinces. For example, in the HB, JS and XJ provinces, both R² and d_r are less than 0.5 (Figure 5b,e).

We also calculate the estimated and statistical mean winter wheat yield from 2016 to 2018 in each province, municipality and county, and the spatial distribution of estimated and statistical mean winter wheat yield has a good consistency over all three scales (Figure 7). Overall, the estimated and statistical mean winter wheat yield are higher in the SD, HN and HB than in other provinces, which also shows at municipal and county scales (Figure 7). Otherwise, in HuB, SC and GS, the mean yield has lower values at all three scales. Moreover, the difference between estimated and statistical mean yield is only obvious at the county scale in XJ (Figure 7(c1,c2)). In addition, at the province scale, the model explains 81% of the variation in yield with a slope of 0.99, and the RMSE, MAPE and d_r are 23.2 g C m⁻² yr⁻¹, 8.8% and 0.8, respectively, which reveals that the estimated yield is consistent with the statistical yield in the corresponding province (Figure 8a). The estimated total production of all investigated 11 provinces is also in good agreement with the statistical value from 2016 to 2018 (Figure 8b). This is especially true in 2016 and 2017, in which the deviations of the estimated and statistical total production are only approximately 3% and 2%, respectively. All of these results demonstrate the reliability of our method in estimating the winter wheat yield.

3.3. Harvest Index Distribution

Using the statistical yield data from 2016 to 2018, we inverted the harvest index with any two years of statistical data and estimated the GPP, and validated the model performance in estimating the yield of the remaining year. We repeated this procedure three times and each year was used as an independent validation. The mean HI values of each province, municipality and county from 2016 to 2018 are shown in Figure 9, which displays a large spatial heterogeneity. Specifically, at the province scale, the HI values of most provinces are in the range of 0.2–0.3 and 0.3–0.4, both totals account for 72.72%, followed by 0.4–0.5 and 0.5–0.6, with a proportion of 18.18% and 9.09%, respectively (Figure 9a). The HI values of most municipalities also lie in the range of 0.2–0.3 and 0.3–0.4, which corresponds to a proportion of 64.38%, followed by 0.4–0.5 and 0.5–0.6, with a proportion of 23.97% and 10.27%, respectively (Figure 9b). In addition, there exists a very low percentage (<0.7%) of municipalities with HI values less than 0.2 or greater than 0.6 (Figure 9b). The distribution of HI at the county scale is quite similar to that at the municipal scale (Figure 9c). The notable proportion of HI is in the range of 0.3–0.4 and 0.4–0.5, which accounts for 58.8%, followed by 0.2–0.3 and 0.5–0.6, with a proportion of 20.6% and 12.81%, respectively (Figure 9c). The HI values less than 0.2 or greater than 0.6 only account for 4.65% and 3.14%, respectively (Figure 9c).

4. Discussion

As one of the largest producers of wheat, China plays an important role in determining the global supplies of winter wheat and retaining national and global food security [65,66]. However, there are few methods for estimating yield and quantifying production of winter wheat at the national scale. This study developed a new method based on a satellite-based light use efficiency model (i.e., EC-LUE), and demonstrated its credible performance for simulating the spatial variations of winter wheat yield over 11 provinces, which accounted for more than 98% of the production of winter wheat in China. Furthermore, our estimates are performed at a 30-m spatial resolution based on Landsat and Sentinel images and provide an opportunity to investigate the spatial variations at a fine resolution.

The EC-LUE model has been successfully applied to simulate crop GPP and yield at 36 EC crop sites [38]. Moreover, a recent study showed that the EC-LUE model had the ability to simulate the spatial and temporal variabilities of winter wheat yield at a state scale in US [41]. Here, our study further confirmed the feasibility of the model in estimating winter wheat yield and production at a national scale. Thus, it is promising to apply the EC-LUE model to predict the spatial and temporal variabilities of other crop yield on a large scale.

This study used a specific harvest index for each province, municipality and county to convert the estimated GPP by EC-LUE to crop yield. The validation shows that our method can simulate the spatial variations of yield and production. The harvest index has been long established as an essential parameter for estimating crop yield [60,67,68]. Prevailing crop growth models, for example AQUACROP [69], also use the harvest index to estimate crop yield. However, the harvest index is affected by crop varieties, growing conditions and tillage practices [70,71,72,73], which poses a great challenge to the identification of the harvest index on a large spatial scale [74,75]. Numerous studies have made great efforts to estimate the harvest index at the regional scale [71,74,76]. Previous studies also found a strong correlation between the harvest index and the fraction of crop transpiration after anthesis [77,78,79,80]. Subsequently, several efforts have been made to estimate the harvest index by calculating the fraction of crop transpiration after anthesis [75]. However, there are large uncertainties in crop transpiration that further result in estimation errors of the harvest index. Another method for estimating the harvest index is to use parameter inversion by calculating the ratio of statistical crop yield to predicted crop biomass as the mean of the harvest index at the regional scale [81]. A recent study used winter wheat yield data and estimated GPP to invert the harvest index in Kansas, USA [41]. This study also used the method to invert the harvest index of winter wheat at each province, municipality and county and used it to estimate the corresponding yield and production. The inverted harvest index mainly varies from 0.2 to 0.6 at all administrative units, which is comparable to field measurements. For example, Hay [61] found that the harvest index of modern wheat ranged from 0.3 to 0.6; Shearman et al. [82] recorded a harvest index of 0.48 to 0.5 for the best winter wheat in the UK. For dryland wheat in Australia, the harvest index varied from 0.08 to 0.56 [83], and for the semiarid areas in China, the harvest index varied from 0.28 to 0.56 [84]. In addition, He et al. [40] summarized that the range of HI values for winter wheat was from 0.33 to 0.53; these values were derived from converting biomass, gross and net primary product to crop yield. In addition, in this study, the inverted harvest index less than 0.2 or greater than 0.6 occurs in some administrative units. The regions with a harvest index less than 0.2 are mainly concentrated in HuB, SAX, SC and the southern region of GS (Figure 9), which is consistent with the spatial distribution of statistical yield (Figure 7). The statistical winter wheat yield in these areas ranges from 2×10³ kg/ha to 4×10³ kg/ha, which is significantly lower than the yield of 6×10³ kg/ha to 8×10³ kg/ha in the main production provinces (HN and SD). The areas with a harvest index greater than 0.6 are mainly distributed in XJ and the northern region of GS.

Several potential causes may result in uncertainties in the GPP and yield estimation. First, we used Landsat7, Landsat8 and Sentinel 2 to synthesize NDVI data. Although these products have been processed for atmospheric correction, the differences of vegetation index resulting from different sensors would affect the quality of NDVI data, and alleviating the influence from this difference remains a challenge [40,57]. Then, the quantity of cloud-free satellite data largely impacts the estimation accuracy of GPP throughout the growing season of winter wheat. This study used the composited images from the Sentinel and Landsat datasets in order to achieve a large fraction of effective satellite data. However, there are large differences in the available satellite images among investigated periods and various provinces and there is a very low fraction of cloud-free images in 2016 and in South China [85]. A recent study has highlighted that the low quantity of effective satellite data largely reduces the accuracy of identifying the planting areas of crops [86]. Here, the model accuracy for simulating yield at the county scale is worse than that at the municipality scale. One potential cause is that the accuracy of identification maps of winter wheat at the county scale is lower when compared to the municipality scale [57].

Although field management, such as irrigation and fertilization, may impact crop photosynthetic capacity, satellite-based vegetation index can usefully indicate the impacts of field management [87]. A previous study demonstrated that the satellite-based LUE model had the capacity to simulate GPP across different field management practices and various crops [38]. However, there are still some uncertainties around several parameters for converting GPP estimations to crop yield, including the ratio of root to aboveground biomass (RS), autotrophic respiration (AR) and moisture content (MC). We inverted the harvest index for counties, municipalities and provinces, but kept the other three parameters constant (i.e., RS = 0.2, AR = 0.66, MC = 0.11). These three parameters may vary with environmental conditions [88]. For example, the ratio of root to aboveground biomass varies with soil moisture [89], soil nutrient [90] and crop variety [91]. Moreover, the regulation of assimilated carbon allocation among crop organs is still poorly understood [92,93,94,95]. Many efforts have been made to investigate how the assimilates are partitioned among plant organs [96,97,98]. However, there is still no validated theory available, and the ratio of root to aboveground biomass is one of the weakest features of current crop growth models [99]. The same issue exists in terms of autotrophic respiration fraction, which usually varies with temperature and growth stages [100,101,102]. However, there are no satisfying methods to simulate the spatial and temporal variations of the fraction of autotrophic respiration [103]. In the future, a higher availability of effective satellite data, a better accuracy of crop planting areas identification and developing methods for describing carbon allocation and autotrophic respiration will help to improve the crop yield and production estimation.

5. Conclusions

This study used the EC-LUE model to estimate vegetation gross primary production (GPP) at a 30 m spatial resolution and convert it to yield by combining with the harvest index and then quantifying the annual winter wheat production at the national scale. To examine the model performance, we first assessed the accuracy of the EC-LUE model to simulate temporal and spatial variations of GPP based on eddy covariance measurements at three winter wheat fields, and then compared the estimated yield with field-surveyed yield and statistical yield at the province, municipal and county scales. The results demonstrate that our method can simulate more than 80% of the spatial variation of GPP at flux sites and more than 50% of the yield spatial variation both in surveyed fields and most administrative units. This is especially true when the model is combined with the region-specific harvest index, and makes it possible to effectively estimate winter wheat yield from county to nation scale. Therefore, our study provides a robust method for converting satellite-based GPP estimates to crop yield on a large spatial scale and would be beneficial for national winter wheat production monitoring.