Winter Wheat Production Estimation Based on Environmental Stress Factors from Satellite Observations

The rapid and accurate estimation of wheat production at a regional scale is crucial for national food security and sustainable agricultural development. This study developed a new gross primary productivity (GPP) estimation model (denoted as the [ACPM]), based on the effects of light, heat, soil moisture, and nitrogen content (N) on the light-use efficiency of winter wheat. The ACPM model used the quantic additivity of the environmental factors to improve the minimum form or multiple multiplication form in the previous model and thus characterized the joint effects of heat, soil moisture, and N on crop photosynthesis performance. The key parameters (i.e., light) were determined from the photosynthetically active radiation product of the Himawari-8 sensor and the fraction of photosynthetically active radiation product of Moderate Resolution Imaging Spectroradiometer (MODIS). The heat was determined from the land temperature products of MODIS. The soil moisture was obtained from the inversion using a visible and shortwave infrared drought index (VSDI), whereas the N stress of winter wheat was detected using the newly developed modified ratio vegetation index (MRVI), which could accurately obtain the spatiotemporal distribution of the leaf chlorophyll content of winter wheat. The ACPM and two other previous models (named the GPP1 and GPP2 models) were applied on the Himawari-8 and MODIS images in Hengshui City. The evaluation results, based on the ground measurement, indicated that the ACPM models exhibited the best estimate of dry aboveground biomass (DAM) and the wheat yield in Hengshui City, with errors of <10% and <12% for the DAM and yield, respectively. Considering the easy operation of the ACPM model and the accessibility of the corresponding satellite images, the Agriculture Crop Photosynthesis Model (ACPM) can be expected to provide information on the winter wheat shortfalls and surplus ahead of the availability of official statistical data.


Study Area
The study area was at Hengshui City (37 • 03 N-38 • 23 N, 115 • 10 E-116 • 34 E) of Hebei Province and included 11 counties, such as the Shenzhou, Zaoqiang, Gucheng, Wuyi, Anping, Raoyang, Wuqiang, Fucheng, Jingxian, Jizhou, and Hengshui counties (Figure 1). The Hengshui City was a part of the North China Plain, which was an alluvial plain of the Yellow River, Huaihe River, and Haihe River [40]. In this region, winter wheat was the dominant crop and the main crop planting patterns were the winter wheat and summer maize rotations. Every year, the winter wheat was sowed in early October and harvested in early June. The climate was semi-humid and semi-arid with an abundant solar heat resource. The mean annual precipitation was about 510 mm, and >70% fell in the summer. Spring was a dry season. The annual mean temperature was 12.7 • C and groundwater depth was below 10 m. The soil was loam fluvo aquic soil. The annual irrigation was about two or three times and the N fertilizer was about 240 kg N/ha for the winter wheat cultivation. The main factor affecting the winter wheat growth was a water deficit with a heavy exploitation of the deep groundwater [41]. In addition, the government conducted a strategy for reducing the fertilizer rate in recent years. Taking into consideration the influence of hydrothermal and fertilizer factors on the crop growth simultaneously, the production-estimation was practical and essential for the winter wheat cultivation in this region.

Study Area
The study area was at Hengshui City (37°03′N-38°23′N, 115°10′E-116°34′E) of Hebei Province and included 11 counties, such as the Shenzhou, Zaoqiang, Gucheng, Wuyi, Anping, Raoyang, Wuqiang, Fucheng, Jingxian, Jizhou, and Hengshui counties (Figure 1). The Hengshui City was a part of the North China Plain, which was an alluvial plain of the Yellow River, Huaihe River, and Haihe River [40]. In this region, winter wheat was the dominant crop and the main crop planting patterns were the winter wheat and summer maize rotations. Every year, the winter wheat was sowed in early October and harvested in early June. The climate was semi-humid and semi-arid with an abundant solar heat resource. The mean annual precipitation was about 510 mm, and >70% fell in the summer. Spring was a dry season. The annual mean temperature was 12.7 °C and groundwater depth was below 10 m. The soil was loam fluvo aquic soil. The annual irrigation was about two or three times and the N fertilizer was about 240 kg N/ha for the winter wheat cultivation. The main factor affecting the winter wheat growth was a water deficit with a heavy exploitation of the deep groundwater [41]. In addition, the government conducted a strategy for reducing the fertilizer rate in recent years. Taking into consideration the influence of hydrothermal and fertilizer factors on the crop growth simultaneously, the production-estimation was practical and essential for the winter wheat cultivation in this region.

Remote Sensing Data
Photosynthetically active radiation (PAR) data from the Himawari-8 satellite, reflectance, land surface temperature (LST), and the FPAR datasets from the MODIS sensors were employed during the winter wheat growth period, from March to June 2017. The Himawari-8 spacecraft was launched on 7 October 2014 from the Tanegashima Space Center (TNSC) of JAXA, Japan. The Himawari standard data included the visible to infrared radiances (Band 1 to 16) of the three regions of the Full

Remote Sensing Data
Photosynthetically active radiation (PAR) data from the Himawari-8 satellite, reflectance, land surface temperature (LST), and the FPAR datasets from the MODIS sensors were employed during the winter wheat growth period, from March to June 2017. The Himawari-8 spacecraft was launched Remote Sens. 2018, 10, 962 5 of 25 on 7 October 2014 from the Tanegashima Space Center (TNSC) of JAXA, Japan. The Himawari standard data included the visible to infrared radiances (Band 1 to 16) of the three regions of the Full Disk (global), Japan Area, and Target Area. The Himawari-8 data were released to research communities for creating the products. The products included short wave radiation, PAR, chlorophyll-a, sea surface temperature, aerosol optical thickness, and wild fire. All of the data were in the NetCDF format. With 5 km of spatial resolution and 30 min of temporal resolution, the PAR data were of level 3 and were downloaded with a JAXA Himawari Monitor (http://www.eorc.jaxa.jp/ptree/userguide.html). The downloaded PAR data were geometrically corrected using the ENVI (The Environment for Visualizing Images) software. Then, the PAR data were integrated into 8-day composite products and resampled to 1 km.
The MODIS images provided a low spatial resolution but a high temporal resolution and wide coverage areas. This study used an 8-day composite time-series dataset of the MODIS surface reflectance (MOD09A1), LST (MOD11A2), and FPAR (MCD15A2H), from 20 March to 31 May 2017. The MODIS dataset were downloaded from EARTHDATA (https://search.earthdata.nasa.gov) with a tile of h27v05. The MOD11A2 MODIS images were of a level 3 product with 1 km pixel size. The MOD09A1 MODIS images were of level 3 product with 500 m pixel size. The MCD15A2H MODIS images were of a level 4 product with a 500 m pixel size. MOD09A1 and MCD15A2H data were resampled to 1 km. Atmospheric, geometric, and radiometric corrections were not required for the MODIS level 3 and 4 data, facilitating the user to popularize and spread the data. However, all of the MODIS images were re-projected from Sinusoidal to the World Geodetic System 1984 projection.
The study used Sentinel-2 time-series data to identify the winter wheat area in Hengshui City. Eight images with 10 m of spatial resolution were captured on 29 March and 7 June 2017, corresponding to the Tillering and Harvest stages of the winter wheat. Atmospheric correction were then performed on the Sentinel-2 images.
The one Worldview2 image that was acquired on 13 March 2017 was used to verify the accuracy of Sentinel-2 classification. The radiometric correction and atmospheric correction were performed on the Worldview2 image. Using a false color composite, the winter wheat distribution map covering part of Shenzhou County was obtained.

Field Data
The field experiments ( Figure 1 and Table 1 Figure 1a). The climate condition and soil texture of exp1 were similar to those of exp2. The following seven types of in situ data were measured: soil moisture (θ), leaf water content, spectral reflectance, chlorophyll content (C ab ), total nitrogen content (TN), dry aboveground biomass (DAM), and yield. In exp1, ten treatments were conducted with five N levels (0, 70, 140, 210 and 280 kg N/ha, which referred to N0, N1, N2, N3 and N4, respectively) and two irrigation water levels (60% and 80% of the field water capacity) from 27 March to 25 May 2011. Ten treatments were conducted in the same field of 0.33 ha. Every treatment had three replications with random arrangements. There were a total 30 plots in exp1 and every plot area was 5 m × 20 m. The measurements were performed in every plot about once a week in exp1. For each experimental sampling plot of exp2, five smaller sub-plots ( Figure 1c) were deployed along diagonal lines, and measurements were performed thrice in 2017, corresponding to the tillering, anthesis, and senescence stages of the winter wheat. A total of 22 plots with 110 sub-plots were measured in 11 counties of Hengshui City (Figure 1b,c). The θ was determined by the gravimetric method. The leaf water content and DAM were measured via the oven drying method, and the C ab and TN in exp1 were determined by spectrophotometry and the micro-Kjeldahl method, and the C ab in exp2 was determined by SPAD-502 Chlorophyll meter (Konlca Minolta Optics, Inc., Williams Drive Ramsey, NJ, USA). The canopy spectral reflectance was measured by an ASD FieldSpec Pro spectrometer (Analytical Spectral Devices, Boulder, CO, USA). Based on the spectral response function of the MODIS, the measured field spectral reflectance was then converted to four multispectral bands corresponding to the blue, green, red, and near-infrared bands (NIR) of the MODIS. The measured θ, C ab , and DAM data in the exp2 plots were averaged by five measurements from the sub-plots. Moreover, the yield in exp2 was averaged by five productions, which were manually collected at a 1 m × 1 m area of the sub-plots at the harvest.

Analyses of Field-and Satellite-Level Data
The field canopy spectrum of the winter wheat was collected in every plot of exp1. The MODIS spectral response function was obtained in the ENVI software. Based on the spectral response function of the MODIS, the field spectral reflectance was converted to the reflectance of the MODIS blue, green, red, and NIR bands. The converted canopy reflectance of the MODIS were analyzed with the measured C ab and TN of exp1, together.
The measured C ab , DAM, and yield of the five subplots in exp2, were averaged for the ground value in the corresponding pixels of the satellite images during the experimental period.
The satellite images with a different resolution and projection coordinates system were obtained in Hengshui City in 2017. All of the images were re-projected to the World Geodetic System in the 1984 projection. The resolution of the PAR data from the Himawari-8 satellite were resized from 5 km to 1 km, using the nearest neighbor algorithm. The resolution of the MODIS surface reflectance (MOD09A1) data were resized from 500 m to 1 km, using the nearest neighbor algorithm as well. The resolution of the MODIS FPAR (MCD15A2H) data were also resized from 500 m to 1 km, using the nearest neighbor algorithm. The resolution of the MODIS LST (MOD11A2) data was 1 km. The satellite images with the resolution of 1 km were prepared for the following inversion.

Winter Wheat Planting Area Dataset
The identification of the winter wheat area was important and necessary for the model development. The extraction accuracy of the winter wheat directly affected the accuracy of the Remote Sens. 2018, 10, 962 7 of 25 yield estimation [42]. For accurately extracting the winter wheat area in Hengshui City, the study used specific phenological characteristics of the winter wheat and Sentinel-2 time-series data with a high spatial resolution of 10 m.
It was well-known that different crops possessed different phenological characteristics in the same region. In Hengshui City, winter wheat was usually sown in October of the former year, turned green in March as the only crop growing in the field, and was then harvested in June when other crops were all green. The area of wheat, evergreen forest, and grassland was extracted based on the higher NDVI compared with the low-threshold (A) at the tillering stage. The planting area of the winter wheat was then determined based on the lower NDVI compared with the up-threshold (B) at the harvest, when the NDVI of the forest and grassland maintained higher values. On the basis of the phenological difference between the winter wheat and the other crops or vegetation in this area, we built a dual threshold method for winter wheat pixel extraction, Equation (1). The classification principle was as follows: where NDVI 1 is the NDVI value at the tillering stage, NDVI 2 is the NDVI value at harvest, A is the low-threshold at the tillering stage and set to 0.6, and B is the up-threshold at harvest and is set to 0.3. Using Sentinel-2 data that were atmospherically corrected, four images from the same day were mosaicked together to cover the whole study area, and the NDVI was calculated. Using the Sentinel-2 NDVI data, the dual threshold classification method was processed to produce a winter wheat map for Hengshui City at the 10 m resolution. The output produced a winter wheat map for Hengshui City with spatial resolution of 10 m in 2017. Using the winter wheat distribution map from the Worldview2 image, the extraction accuracy of the winter wheat from the Sentinel-2 images was verified. An error matrix for the extraction results was created and is shown in Table 2. We randomly took a total 200 points of the winter wheat field and non-wheat field in the winter wheat map from the Sentinel-2 images. Every point in the winter wheat map from the Sentinel-2 images was the size of 10 m × 10 m and 5 × 5 pixels, which corresponded to the winter wheat map from the Worldview2 image. If the proportion of the winter wheat area in the 5 × 5 pixels was higher than 0.6 in the Worldview2 image, the corresponding pixel in the winter wheat map from the Sentinel-2 images was a point of the winter wheat field. As shown in Table 2, the overall accuracy of the classification was 93.5%.
For scale consistency, the binary mask was resized to produce a final area mask with 1 km resolution by computing the percent of the wheat value (Pct) from the Sentinel-2 binary mask (10 m resolution) map within each 1 km pixel footprint ( Figure 2) [3]. When the Pct was >80%, the pixel was pure for the winter wheat. When the Pct was distributed from 0.5 to 0.8, the pixel was a mixed-pixel. If the Pct was <0.5, the pixels were neglected for the verification of the field data below. Six pure pixels and eleven mixed pixels were at the experimental sampling plots of exp2.

The Schema of Research Procedure
The schema of this research procedure was illustrated in Figure 3. Four steps were involved in this procedure. First of all, for detecting the wheat canopy Cab, we developed the new modified ratio vegetation index (MRVI), based on the field canopy reflectance, Cab, and TN data in exp1, at Yucheng Station. Secondly, the key parameters in the GPP estimation models were inversed, for example, MRVI; ScaledWDRVI; ScaledLST; ScaledVSDI (visible and shortwave infrared drought index), Equations (5) to (7); and the green normalized difference vegetation index (GNDVI). Thirdly, the new GPP estimation model (ACPM) was developed for the winter wheat in Hengshui City. Fourthly, the yield and DAM were simulated based on the ACPM, GPP1, and GPP2 models and the accuracies of ACPM for DAM and yield estimation were verified.

The Schema of Research Procedure
The schema of this research procedure was illustrated in Figure 3. Four steps were involved in this procedure. First of all, for detecting the wheat canopy C ab , we developed the new modified ratio vegetation index (MRVI), based on the field canopy reflectance, C ab , and TN data in exp1, at Yucheng Station. Secondly, the key parameters in the GPP estimation models were inversed, for example, MRVI; ScaledWDRVI; ScaledLST; ScaledVSDI (visible and shortwave infrared drought index), Equations (5) to (7); and the green normalized difference vegetation index (GNDVI). Thirdly, the new GPP estimation model (ACPM) was developed for the winter wheat in Hengshui City. Fourthly, the yield and DAM were simulated based on the ACPM, GPP1, and GPP2 models and the accuracies of ACPM for DAM and yield estimation were verified.

The Schema of Research Procedure
The schema of this research procedure was illustrated in Figure 3. Four steps were involved in this procedure. First of all, for detecting the wheat canopy Cab, we developed the new modified ratio vegetation index (MRVI), based on the field canopy reflectance, Cab, and TN data in exp1, at Yucheng Station. Secondly, the key parameters in the GPP estimation models were inversed, for example, MRVI; ScaledWDRVI; ScaledLST; ScaledVSDI (visible and shortwave infrared drought index), Equations (5) to (7); and the green normalized difference vegetation index (GNDVI). Thirdly, the new GPP estimation model (ACPM) was developed for the winter wheat in Hengshui City. Fourthly, the yield and DAM were simulated based on the ACPM, GPP1, and GPP2 models and the accuracies of ACPM for DAM and yield estimation were verified.

The Modified Ratio Vegetation Index (MRVI) for Detecting Wheat Canopy C ab
It was reported that the ratio of NIR to green (RVI2) was linearly related to the total leaf N accumulation and that it was useful for monitoring the rice N status [43]. However, the higher deviation was found using the RVI2 for monitoring the winter wheat C ab in exp1 at Yucheng Station [44]. Thus, the study involved the addition of blue band information to modify the ratio index (RVI2), to improve the inversion accuracy using multispectral images, for example, MODIS. The new modified ratio vegetation index (MRVI) was defined as follows: where ρ NIR is the reflectance at NIR, ρ green is the reflectance at green band, ρ blue is the reflectance at blue band, α is the empirical coefficient, which controls the dispersion of MRVI values and lets the value of MRVI change between 0 and 1, was set to 35 in the study.

GPP Estimation Model Based on LUE and Environmental Factors
Considering the effects of the environmental factors on the GPP accumulation, Zhang N (2013) built two GPP estimation models based on LUE [44]. Like other LUE models, for example, VPM, CASA, MODIS-GPP, and EC-LUE, the models were based on the one minimum or coupled effect of the environmental factors on crop growth, and were designed as follows: where LUE max is the maximum light use efficiency of the crop, 1.02-3.71 gC/MJ for wheat [45]; ScaledWDRVI is the influence of the sunshine; WDRVI is the wide dynamic range vegetation index [46]; ScaledLST is the influence of temperature; LST is the land temperature product from MODIS (MOD11A2); ScaledVSDI is the influence of soil moisture; VSDI is the visible and shortwave infrared drought index [44]; and GNDVI is the influence of N, which is defined as the green normalized difference vegetation index [47]. Solely considering the effect of the worst effect of one factor or the coupled effect of factors was not practical for field crop growth, because light, temperature, soil moisture, and N affected the crop growth simultaneously. The crop growth was the joint result of external environmental factors, and heat, soil moisture, and N played complementary roles in crop growth. A sufficient hydrothermal environment could compensate the yield-reducing effect that was caused by the N deficit. Similarly, an increase in the N rate could improve the yield-reducing effect that was caused by the water deficit.
The GPP accumulation was mainly controlled by photosynthesis, which was determined by the PAR, FPAR, and LUE. The LUE was determined by the LUE max and adjustment coefficients of the N deficit, water droughts, and temperature. Thus, based on the LUE theory, the study modified the aforementioned model and designed it as follows: where ACPM is the Agriculture Crop Photosynthesis Model and FPAR is the product from MODIS (MCD15A2H).
Using the three above-mentioned models, we estimated the GPP for the winter wheat in Hengshui City in 2017. The DAM and yield were then estimated based on the GPP accumulation across the whole growth period of the winter wheat.

DAM and Yield Estimation for Winter Wheat
Neglecting the differences of species, we assumed the same transformation process from photosynthate to DAM or yield in all of winter wheat fields. Considering the available relationship between GPP and DAM and yield for wheat [48,49], we determined the winter wheat DAM and yield as follows: where CUE is the carbon use efficiency and is set to 0.5 [50], HI is the harvest index and is set to 0.45 [49], RSR is the belowground to aboveground biomass ratio and is set to 0.2 [51], θ is the grain moisture at harvest and is set to 11% [49], and CR is the carbon content and is set to 0.45 [49]. Based on the three aforementioned GPP estimation models, the winter wheat DAM and yield were evaluated for Hengshui City in 2017. The optimized model would achieve better forecasting accuracy for the DAM and yield estimation, prior to the government obtaining the statistical field yield data.

Verificaiton of Model and New Index
Based on the assumption of the consistent relationship between the N accumulation and the canopy reflectance and ratio index across the whole growth period, the VIs (listed in Table 3) were collected from the previous literature [43,47,[52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69]. Statistical analyses were performed using SPSS 13.0. Linear correlation analyses between the C ab and VIs were conducted individually to evaluate the feasibility of each vegetation index. The precision of the VIs was evaluated using the square of the correlation coefficient of the linear regression (R 2 ) and the root mean square error (RMSE). Table 3. The relationship between vegetation indexes (VIs) and the leaf C ab of the winter wheat in exp1 at the Yucheng Comprehensive Experimental Station.

Indexes
Expressions  To validate the accuracy of the GPP models, the DAM and yield of the winter wheat was simulated in Hengshui City in 2017. The ground data of the DAM and yield in the corresponding pixels were used to validate the estimated results. Based on the linear regression analysis with the intercept of zero, the performance of the model was assessed using the correlation coefficient of the linear regression (r), RMSE, and the error. The error was calculated as follows: where RMSE is the root mean square error and X is the mean value of the measured data of the DAM and yield of the winter wheat.

Results
The corresponding leaf, TN, and canopy spectrum with different N levels were analyzed during the growth period of the winter wheat in exp1, as well as the accuracy of the MRVI for detecting the leaf C ab were validated in exp1 and exp2, respectively. The GPP estimation models were run for each winter wheat pixel (1 km) in Hengshui City in 2017. The situ yield and DAM, with the estimated values in the corresponding pixel, were compared in exp2 and the accuracy of the models for the yield estimation was validated.

Variation of Leaf Chlorophyll Content (C ab ), Total Nitrogen Content (TN), and Canopy Spectrum at Different N Level
As shown in Figure 4, the winter wheat leaf C ab exhibited a significant linear relationship with the leaf TN distribution, with a high correlation coefficient of 0.8 in exp1. It demonstrated that the variation of the leaf C ab that exposed the winter wheat leaf TN trend during the whole winter wheat growth period.
The variation of the winter wheat canopy reflectance with wavelength at a different date and with different N levels is shown in Figure 5a-d. The variation of leaf C ab with time at different N levels was shown in Figure 5e. As shown in Figure 5, the winter wheat canopy reflectance and leaf C ab all varied with different N rates during the winter wheat growth period, especially at the NIR band. The reflectance at the NIR band obviously increased with the increase of the N level, but the difference of the reflectances at the NIR band gradually decreased with N level >140 kg N/ha. Furthermore, the leaf C ab also increased with the increase of the N level, and the difference of the leaf C ab under different treatments was still obvious at 140-280 kg N/ha.
The variation of the winter wheat canopy reflectance with the leaf C ab is shown in Figure 6. Linear correlation analyses were conducted and the correlation coefficients were significant (p < 0.01), shown in Figure 6a-f. As shown in Figures 5 and 6b, the reflectance at the NIR band was the most sensitive to different N levels and the variation of C ab , with the highest R 2 being 0.49 and lowest RMSE being 9.47 µg/cm 2 . The reflectance at the red and green bands were more sensitive to the variation of C ab , but tended to change little when the C ab was larger than 40 µg/cm 2 (Figure 6a,d). The reflectance at the blue band was lower and tended to change only a little as well, when C ab was larger than 40 µg/cm 2 . For highlighting the difference of the red, blue, and green bands, we used the mathematical transformation. As shown in Figure 6e-h, the linear relationships between ρ green / ρ green − ρ blue 2 , ρ blue / ρ green − ρ blue 2 and the leaf C ab were higher with the R 2 of > 0.5.
By contrast, the ρ red / ρ red − ρ green 2 and ρ green / ρ red − ρ green 2 had several outliers because the lower reflectance in the red band and the lower R 2 of <0.2 were obtained (p > 0.05). Furthermore, ρ blue / ρ green − ρ blue 2 showed a lower dispersion than ρ green / ρ green − ρ blue 2 , with the RMSE of 7.97.
The new vegetation index served the new algorithms of the GPP estimation. The lower dispersion of results, the better. Thus, we used the NIR band and ρ blue / ρ green − ρ blue 2 to develop the new vegetation index (MRVI) for monitoring the C ab of winter wheat.

Validation of the Modified Ratio Vegetation Index (MRVI)
The existing VIs for monitoring the crop N content or C ab were collected from previous literatures (Table 3) [43,47,[52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69]. Based on the MODIS reflectance at the blue, green, red, and NIR bands, which were transformed from the field spectral reflectance in exp1, the VIs were calculated and are listed in Table 3. The relationships between the VIs and leaf C ab in exp1 were discussed using the linear regression analysis. The VI was the x-axis and the C ab was the y-axis. The accuracy of VI was evaluated using a correlation coefficient (R 2 ) and the root mean square error (RMSE). If the R 2 was higher and RMSE was lower, the VI was recognized to be more accurate for detecting the C ab of winter wheat. As shown in Table 3, the results showed that the MRVI obtained the best estimation accuracy with an increased correlation coefficient of 0.62 and it decreased the RMSE of 8.34 µg/cm 2 . The improved results were followed by RVI2, a newly built index (ρ NIR /ρ blue ); MSR; and SIPI, with correlation coefficients of 0.59, 0.54, 0.52 and 0.52, respectively, corresponding to RMSEs of 8.66, 9.08, 9.27, and 9.35 µg/cm 2 , respectively. Except for normalized difference vegetation index (green-blue) (NDVIg-b) and plant senescence reflectance index (PSRI), the correlation coefficients of C ab and VIs were significant (p < 0.05). The MRVI improved the estimation accuracy of the leaf C ab of the winter wheat in comparison with other indexes.   The relationship between the canopy reflectance and leaf C ab of winter wheat in exp1 at Yucheng Station, as follows: (a) the reflectance at the red band and leaf C ab ; (b) the reflectance at the near-infrared (NIR) band and leaf C ab ; (c) the reflectance at the blue band and leaf C ab ; (d) the reflectance at the green band and leaf C ab ; (e) ρ blue / ρ green − ρ blue 2 and the leaf C ab ; (f) ρ green / ρ green − ρ blue 2 and the leaf C ab ; (g) ρ red / ρ red − ρ green 2 and the leaf C ab ; and (h) ρ green / ρ red − ρ green 2 and the leaf C ab .
The verified five indexes with the higher R 2 and lower RMSE in exp1, that is, MRVI, RVI2, ρ NIR /ρ blue , MSR, and SIPI, were calculated for each winter wheat pixel (1 km) in Hengshui City in 2017. The relationships between the situ leaf C ab and the calculated indexes of the corresponding pixel were analyzed using the linear regression in exp2 and the accuracy of the index was validated using the R 2 . There were six pure pixels and eleven mixed pixels at the experimental plots of exp2. The analyzed results of the leaf C ab and VIs for the pure pixels and all of the pixels of exp2, are shown in Figures 7 and 8, respectively. Figure 7 shows the relationships between the leaf C ab and the five indexes (MRVI, RVI2, ρ NIR /ρ blue , MSR, and SIPI) in the pure winter wheat pixels of exp2, at Hengshui City. Figure 8 shows the relationships between the leaf C ab and the five indexes in all of the winter wheat pixels of exp2 (pure pixels and mixed pixels). As shown in Figures 7 and 8, during the winter wheat growth period of Hengshui City, the top leaf C ab was well estimated by MRVI, with higher correlation coefficients, which were both significant (p < 0.01). The MRVI obtained the highest correlation coefficient of 0.73 in the pure pixels (Figure 7a), whereas the others were about 0.30-0.39 (Figure 7b-e). Furthermore, as shown in Figure 8a, the MRVI also obtained the highest correlation coefficient of 0.59 in all of the pixels (including the pure and mixed pixels), whereas the others were about 0.23-0.31, which are shown in Figure 8b-e. There were outliers in the mixed pixels, and the reason for that was mainly that the reflectance of the MODIS pixel was susceptible to the mixed-pixel effect with different land cover types. Although the mixed-pixel reduced the estimation accuracy of VIs, the MRVI obtained the better estimation results of the leaf C ab in Hengshui City. It was thus acceptable that the MRVI represented the effect of N on the winter wheat growth, instead of the GNDVI in the GPP1 and GPP2 models.
There were six pure pixels and eleven mixed pixels at the experimental plots of exp2. The analyzed results of the leaf Cab and VIs for the pure pixels and all of the pixels of exp2, are shown in Figures 7  and 8, respectively. Figure 7 shows the relationships between the leaf Cab and the five indexes (MRVI, RVI2, ρ NIR /ρ blue , MSR, and SIPI) in the pure winter wheat pixels of exp2, at Hengshui City. Figure 8 shows the relationships between the leaf Cab and the five indexes in all of the winter wheat pixels of exp2 (pure pixels and mixed pixels). As shown in Figures 7 and 8, during the winter wheat growth period of Hengshui City, the top leaf Cab was well estimated by MRVI, with higher correlation coefficients, which were both significant (p < 0.01). The MRVI obtained the highest correlation coefficient of 0.73 in the pure pixels (Figure 7a), whereas the others were about 0.30-0.39 (Figure 7be). Furthermore, as shown in Figure 8a, the MRVI also obtained the highest correlation coefficient of 0.59 in all of the pixels (including the pure and mixed pixels), whereas the others were about 0.23-0.31, which are shown in Figure 8b-e. There were outliers in the mixed pixels, and the reason for that was mainly that the reflectance of the MODIS pixel was susceptible to the mixed-pixel effect with different land cover types. Although the mixed-pixel reduced the estimation accuracy of VIs, the MRVI obtained the better estimation results of the leaf Cab in Hengshui City. It was thus acceptable that the MRVI represented the effect of N on the winter wheat growth, instead of the GNDVI in the GPP1 and GPP2 models.

Evaluation of DAM of Winter Wheat at Field Scale
The DAM in Hengshui City was evaluated based on the three GPP estimation models at the mature period of the winter wheat in 2017 (Figures 9 and 10). The difference in the winter wheat DAM was only due to the different GPP estimation models. The results showed that the evaluated DAM based on GPP1 and GPP2 estimation was <6 t/ha (Figure 9a,b). The evaluated DAM, based on the GPP1 was the lowest and the evaluated DAM, based on GPP2 was also much lower than the actual DAM of the winter wheat (Figure 10a-d). The main reason was that the serious

Evaluation of DAM of Winter Wheat at Field Scale
The DAM in Hengshui City was evaluated based on the three GPP estimation models at the mature period of the winter wheat in 2017 (Figures 9 and 10). The difference in the winter wheat DAM was only due to the different GPP estimation models. The results showed that the evaluated DAM based on GPP1 and GPP2 estimation was <6 t/ha (Figure 9a,b). The evaluated DAM, based on the GPP1 was the lowest and the evaluated DAM, based on GPP2 was also much lower than the

Evaluation of DAM of Winter Wheat at Field Scale
The DAM in Hengshui City was evaluated based on the three GPP estimation models at the mature period of the winter wheat in 2017 (Figures 9 and 10). The difference in the winter wheat DAM was only due to the different GPP estimation models. The results showed that the evaluated DAM based on GPP1 and GPP2 estimation was <6 t/ha (Figure 9a,b). The evaluated DAM, based on the GPP1 was the lowest and the evaluated DAM, based on GPP2 was also much lower than the actual DAM of the winter wheat (Figure 10a-d). The main reason was that the serious underestimation of the GPP2 and GPP1 models resulted in the decreased DAM accumulation. By contrast, the ACPM model obtained improved estimation accuracy with an error of <10% (Figure 10e,f). Based on the results of the ACPM model in Hengshui City, the RMSE of the measured DAM versus the estimated DAM was 0.98 t/ha, which was equivalent to a 7% error in the pure winter wheat pixels (Figure 10e). Setting the intercept to zero, the linear regression slope was 1.0141, and the correlation coefficient was 0.47 in the pure pixels. Furthermore, the RMSE of the measured DAM versus the estimated DAM, by the ACPM model, was 1.33 t/ha, which was equivalent to a 9% error in the pure and mixed winter wheat pixels (Figure 10f). The linear regression slope was 1.0205, and the correlation coefficient was 0.59 in the pure and mixed pixels. The mixed-pixel increased the estimation error, but the ACPM clearly obtained a higher estimation accuracy of the DAM than the GPP1 and GPP2 models.

Evaluation of Winter Wheat Yield at Field Scale
The winter wheat yield was evaluated at Hengshui City in 2017 and the results are shown in Figures 11 and 12. Similarly, the models that were based on the GPP2 and GPP1 estimation significantly underestimated the yield of the winter wheat (Figures 11 and 12). The evaluated yield based on GPP1 was the lowest, with an error of 92% (Figure 12a,b). The evaluated yield based on GPP2 was also much lower than the actual DAM of the winter wheat, with an error of 80% ( Figure  12c,d). By contrast, the model based on ACPM obtained the best accuracy for yield estimation. Based on the ACPM estimation results, the RMSE of the measured yield versus the estimated yield was 0.51 t/ha, which was equivalent to a 8% error in the pure pixels (Figure 12e). Setting the intercept to zero, we observed a linear regression slope of 1.0126 and a regression coefficient of 0.59 in the pure pixels. In addition, the RMSE of the measured yield versus the estimated yield was 0.85 t/ha, which was equivalent to a 12% error in the pure and mixed pixels (Figure 12f). Setting the intercept to zero, we

Evaluation of Winter Wheat Yield at Field Scale
The winter wheat yield was evaluated at Hengshui City in 2017 and the results are shown in Figures 11 and 12. Similarly, the models that were based on the GPP2 and GPP1 estimation significantly underestimated the yield of the winter wheat (Figures 11 and 12). The evaluated yield based on GPP1 was the lowest, with an error of 92% (Figure 12a,b). The evaluated yield based on GPP2 was also much lower than the actual DAM of the winter wheat, with an error of 80% (Figure 12c,d). By contrast, the model based on ACPM obtained the best accuracy for yield estimation. Based on the ACPM estimation results, the RMSE of the measured yield versus the estimated yield was 0.51 t/ha, which was equivalent to a 8% error in the pure pixels (Figure 12e). Setting the intercept to zero, we observed a linear regression slope of 1.0126 and a regression coefficient of 0.59 in the pure pixels. In addition, the RMSE of the measured yield versus the estimated yield was 0.85 t/ha, which was equivalent to a 12% error in the pure and mixed pixels (Figure 12f). Setting the intercept to zero, we obtained a linear regression slope of 1.0401 and a regression coefficient of 0.46 in the pure and mixed pixels. The mixed pixel reduced the accuracy of the yield estimation, but the results demonstrated that the yield estimation model, based on the ACPM model, was available and practical for Hengshui City. Furthermore, as shown in Figure 10, the Shenzhou County was the highest-yield region, while the Zaoqiang and Wuyi Counties were the lowest-yield regions. These results were very consistent with the actual distribution.
on the ACPM estimation results, the RMSE of the measured yield versus the estimated yield was 0.51 t/ha, which was equivalent to a 8% error in the pure pixels (Figure 12e). Setting the intercept to zero, we observed a linear regression slope of 1.0126 and a regression coefficient of 0.59 in the pure pixels. In addition, the RMSE of the measured yield versus the estimated yield was 0.85 t/ha, which was equivalent to a 12% error in the pure and mixed pixels (Figure 12f). Setting the intercept to zero, we obtained a linear regression slope of 1.0401 and a regression coefficient of 0.46 in the pure and mixed pixels. The mixed pixel reduced the accuracy of the yield estimation, but the results demonstrated that the yield estimation model, based on the ACPM model, was available and practical for Hengshui City. Furthermore, as shown in Figure 10, the Shenzhou County was the highest-yield region, while the Zaoqiang and Wuyi Counties were the lowest-yield regions. These results were very consistent with the actual distribution.

Discussion
The identification of winter wheat was important for yield forecasting, especially for the government county-level crop statistics [70]. Remote sensing data with a high spatial resolution, such as 10 or 30 m, was used for the extraction of the crop area. The study extracted the winter wheat area of Hengshui City from Sentinel-2 images with an overall accuracy of 93.5%. Based on the MODIS EVI data and China's Environment Satellite (HJ-1) data, Yao et al. [71] extracted the corn area over a regional scale with the support vector machine algorithm and the overall accuracy was 82.17%. This result was a little lower than our result of classification. Using the support vector machine algorithm as well, Zhang et al. [72] obtained the classification for corn, based on the MODIS EVI data and HJ-1 images over the whole of Northeast China, and the overall accuracy was 79%, which was also a little lower than our result. Our classification result was available for a distribution map of the winter wheat yield.
The higher deviation of environmental factors led to the increased error of yield estimation. The goal of the newly developed MRVI was reducing the RMSE and increasing the accuracy of C ab estimation. The results indicated that leaf C ab could be estimated using MODIS MRVI data with an accuracy of 8.34 µg/cm 2 . Delegido et al. [73] obtained an RMSE of 4.2 µg/cm 2 for the leaf C ab estimation, using the normalized area over the reflectance curve (NAOC) from the hyperspectral data of the airborne CASI data. Vincini et al. [74] reported an RMSE of 4-5 µg/cm 2 for the leaf C ab estimation with different canopy structures, using the different indexes, e.g., TCI/OSAVI and MTCI from the Sentinel-2 data as well. Those results were both better than our result. However, different sensors had different inversion effects. The higher resolution of the airborne and Sentinel-2 images could improve the accuracy. Using the physical model (The REGularized canopy reFLECtance model), Houborg et al. [75] obtained the higher accuracies with the RMSE of 4.4 µg/cm 2 , based on the 1 m resolution aircraft and with an RMSE of 7.1 µg/cm 2 based on the 10 m resolution SPOT-5 imagery. However, this approach was time-consuming. Li et al. [76] also reported an RMSE of 9.35-13.36 µg/cm 2 for the leaf C ab estimation, using hyperspectral indices, which were comparable with our result. Zarco-Tejada et al. [69] also obtained an RMSE of 11.5 µg/cm 2 for the leaf C ab estimation, using the optimized soil-adjusted vegetation index (OSAVI), and the RMSE was slightly higher than our result. Furthermore, in comparison with the existing VIs from the previous literature, listed in Table 3, the MRVI obtained a better accuracy for the leaf C ab estimation in exp1 at the Yucheng Station. For detecting the leaf C ab distribution in Hengshui City in 2017, MRVI showed better results as well. Thus, it was proper that the N stress of the winter wheat was detected using a newly developed MRVI for the GPP estimation.
Improved results of DAM and yield estimation were obtained in Hengshui City in 2017 after utilizing the modified GPP estimation model (ACPM). Compared with the DAM and yield results based on GPP1 and GPP2 models, the results based on the ACPM model were significantly increased. Compared with the field data, the simulated DAM and yield results were consistent with the actual distribution in Hengshui City. The error of DAM was <10%, and the error of yield was <12%. Similarly, based on the GPP estimation in Boreal Ecosystem Productivity Simulator (BEPS), Wang et al. [77] simulated the winter wheat yield with mean relative error of 4.6% in the North China Plain, but this model hardly satisfied the practical production, owing to many input parameters (e.g., weather condition). A biomass relative error of 28% was reported in Southwest France after coupling formosat-2 data with a light-use-efficiency algorithm for the yield estimate [36]. In comparison with this result, the DAM results in the current study showed an improved accuracy to some extent. Simultaneously, other empirical models from the MODIS NDVI or EVI time series data forecasted a winter wheat yield in the USA or Hungary within 10% of the official reported yield, 2-6 weeks prior to harvest or 7-10 days before the corn silking stage, which was comparable with our results [3,6,26]. However, those approaches were hard to apply in other regions. Introducing the crop growth monitoring system into the China Agriculture Remote Sensing Monitoring System, Fei et al. [78] reported a more than 88% accuracy of the yield forecasting for the winter wheat in the North China Plain. This result was comparable with our result of the yield estimation. Son et al. [79] performed the comparative analysis of MODIS EVI and NDVI time series data for the rice yield estimation and the EVI-based model obtained the higher accuracy with a RMSE of 6.9-8.1%, which was a little better than our result. However, the accuracy of the forecasting method for county or state areas would become less reliable for small areas, and few literature focussed on the quantitative estimation of the yield at the field scale using the satellite data [6,80]. In the current study, our results were verified by in situ yield data, and the accuracy was comparable to the aforementioned studies with acceptable estimation yields. Zhao et al. [81] used the difference water index (RDWI) from the land satellite thematic map (TM) and the NIR band from the MODIS images to predict the winter wheat yield in the North China Plain, and the average relative error between the situ yield and the estimated yield was 10.11%, which was comparable with our results. Wang et al. [82] compared the performance of three modeling algorithms, namely, the Random Forest (RF) machine-learning algorithm, the RF algorithm with support vector regression (SVR), and artificial neural network (ANN) machine-learning algorithms, for the wheat-biomass estimation based on the 15 VIs from the environment and disaster reduction small satellites (HJ-CCD) images, the RF model obtained the more accurate estimation. Fu et al. [83] compared the accuracy of the narrow band vegetation indices and red-edge position for the winter wheat biomass estimation and the partial least square regression method, based on the combination of the band depth ratio (BDR), and the optimal NDVI-like obtained the best estimation result with a RMSE of 1.77 t/ha, which was clearly higher than our result. Using the step-wise regression, Liaqat et al. [84] compared the accuracy of different indexes (e.g., SAVI, EVI, NDVI and MSAVI from MODIS images) for the wheat yield-estimation in the irrigated Indus Basin, and the SAVI obtained the best estimation result with a R2 of 0.63, which was slightly higher than our result. However, these models were hard to apply in other regions as a result of a lack of the physical mechanism. Using the PSO algorithm, Jin et al. [85] assimilated the combined EVI × RVI data from the HJ-1A/B + RADARSAT-2 images into the AquaCrop for yield estimation of the winter wheat, and the RMSE of the estimated yield with the measured data was 0.81 t/ha, which was comparable with our result. With an ensemble Kalman filter (EnKF) algorithm, Xie et al. [86] compared the accuracies of five assimilation schemes, which assimilated the LAI and soil moisture at different wheat growth stages into the CERES-wheat for the yield estimation. The results showed that the assimilation of LAI at the jointing and heading-filling stages into CERES-wheat obtained the higher accuracy, with a RMSE of 548.97 kg/ha, which was better than our estimated result [86]. Adding the GLDAS/Noah-derived surface incoming solar radiation as an input, MODIS-based LAI was also assimilated using a sequential update algorithm into the Soil Water Atmosphere Plant model, and the results of the yield-estimation were improved by 14-26%, with the absolute errors of <7%, which was better than our result [87]. Huang et al. [11] assimilated the LAI from the TM and MODIS data into the WOFOST model for the winter wheat yield estimation, and the error for yield-estimation was <8%, which was better than our result as well. However, the assimilation algorithms were time consuming and needed much more inputs than the ACPM model (e.g., meteorological data and soil data). Furthermore, the ACPM model was easy to operate and the corresponding satellite images were accessible. Thus, the model could be expected to provide reliable information on the winter wheat shortfalls and surplus.
The yield estimation was performed at the harvest, whereas the forecasts were performed several weeks before the harvest [17]. However, the wheat GPP accumulation was vulnerable to environmental stress during the period from the anthesis to the ripening stage. Our estimation was conducted at the beginning of the harvest of the winter wheat in Hengshui City, and could be an indication of the winter wheat shortfalls and surplus, prior to the statistical data that were reported by the government. In addition, if the interannual time series PAR data from 2000 to 2014 could be obtained in the future, the effect of the interannual climate change on the yield could be considered, and our method could be modified for a forecasting yield.

Conclusions
The rapid and accurate estimation of wheat production is important to national food security and sustainable agricultural development. In this study, we developed a new GPP estimation model named the Agriculture Crop Photosynthesis Model (ACPM) for the winter wheat productivity estimation, based on the light-use efficiency theory and environmental stress factors from the Himawari-8 and MODIS satellite observations. The main conclusions were cited as follows: 1.
We developed the new modified ratio vegetation index (MRVI), which could reveal the detailed spatiotemporal distribution of the leaf C ab and N status of the winter wheat in Hengshui City.

2.
We newly developed the ACPM for the winter wheat productivity estimation based on light-use efficiency theory and environmental stress factors from the Himawari-8 and MODIS satellite observations. This model described the joint effects of heat, soil moisture, and N on the crop photosynthesis performance. The ACPM model used a quantic additivity of the environmental factors in order to improve the minimum form or multiple multiplication form in the previous models. The light was determined from the MODIS FPAR data and Himawari-8 PAR data. The heat was determined from the MODIS LST data. The soil moisture was obtained from the inversion, using a visible and shortwave infrared drought index (VSDI). The N stress of the winter wheat was detected using MRVI.

3.
Based on the newly developed GPP model (ACPM), the DAM and yield were well estimated within a 10% and 12% error of the situ data in Hengshui City in 2017. Comparing the DAM and yield results based on the ACPM, GPP1, and GPP2 models, the ACPM model improved the underestimation of the DAM and yield results, based on the previous GPP1 and GPP2 models.
This new ACPM model was simple, economical, timely, robust, and easily operational. Considering the easy operation of the ACPM model and the accessibility of the corresponding satellite images, the model can be expected to provide information on the winter wheat shortfalls and surplus ahead of the availability of the official statistical data.