A Method for Monitoring and Forecasting the Heading and Flowering Dates of Winter Wheat Combining Satellite-Derived Green-up Dates and Accumulated Temperature

Abstract

Heading and flowering are two key phenological stages in the growth process of winter wheat. It is of great significance for agricultural management and scientific research to accurately monitor and forecast the heading and flowering dates of winter wheat. However, the monitoring accuracy of existing methods based on remote sensing needs to be improved, and these methods cannot realize forecasting in advance. This study proposed an accumulated temperature method (ATM) for monitoring and forecasting the heading and flowering dates of winter wheat from the perspective of thermal requirements for crop growth. The ATM method consists of three key procedures: (1) extracting the green-up date of winter wheat as the starting point of temperature accumulation with the dynamic threshold method from remotely sensed vegetation index (VI) time-series data, (2) calculating the accumulated temperature and determining the thermal requirements from the green-up date to the heading date or the flowering date based on phenology observation samples, and (3) combining the satellite-derived green-up date, daily temperature data, and thermal requirements to monitor and forecast the heading date and flowering date of winter wheat. When applying the ATM method to winter wheat in the North China Plain during 2017–2019, the root mean square error (RMSE) for the estimated heading date was between 4.76 and 6.13 d and the RMSE for the estimated flowering date was between 5.30 and 6.41 d. By contrast, the RMSE for the heading and flowering dates estimated by the widely used maximum vegetation index method was approximately 10 d. Furthermore, the forecasting accuracy of the ATM method was also high, and the RMSE was approximately 6 d. In summary, the proposed ATM method can be used to accurately monitor and forecast the heading and flowering dates of winter wheat in large spatial scales and it performs better than the existing maximum vegetation index method.

1. Introduction

Heading and flowering are two important phenological stages in the growth process of winter wheat. The heading date is the time node that marks the transition from vegetative growth to reproductive growth [1], and the flowering date is a key period that is closely related to the final yield of winter wheat. During the flowering phase, the growth and development of winter wheat is very sensitive to factors such as water and heat conditions, soil fertility, and field management levels [2,3,4], and is a critical stage where agricultural practices and crop yield estimation research are concerned [5,6]. Therefore, accurately monitoring the heading and flowering dates of winter wheat is of great significance to guide agricultural practices and conduct related scientific research.

Currently, crop phenological information is mainly derived through two methods: ground observation and retrieval from satellite time-series data. Ground observations monitor the growth and development status of field crops through direct human observations and records. The advantages of ground observation are rooted in their simplicity, high precision, and long time span [7]. However, this method is labor intensive and cannot provide spatially continuous phenological information at large scales. To solve these problems, the method of retrieving phenology from satellite time-series data is usually used. Satellite remote sensing technology has been investigated and practiced successfully for retrieving crop phenology based on remotely sensed vegetation index (VI) time-series data [8,9,10]. At present, a variety of methods have been developed to extract crop phenology based on remote sensing time-series data, mainly including the threshold method [11,12], derivative method [13], moving average method [14], and function fitting method [15]. However, most of the existing methods only focus on the start of the season (SOS) and the end of the season (EOS), lacking effective methods to monitor other crucial crop phenology such as the heading date and flowering date.

Among the methods that can monitor the heading and flowering dates of winter wheat, the maximum vegetation index method (denoted VI_max) is widely used. Because the VI_max method is easy and fast, many studies have chosen this method to extract the heading date or the flowering date of crops [1,16,17,18,19,20]. Moreover, in recent years, other ideas or methods have been proposed to extract the heading date or the flowering date of winter wheat, such as the two-step filtering (TSF) method [21,22], methods based on radar backscatter coefficient time-series data [10,23,24], and methods through transforming temporal scales [25,26,27].

Although there are some methods to extract the heading date and the flowering date of winter wheat, there is still room for improvements. Firstly, the extraction accuracy could be improved. Taking the VI_max method as an example, the root mean square error (RMSE) of the heading date or the flowering date extracted by the method is 10–15 d [13,16,19,20]. The results may be acceptable in large-scale ecological applications, but it is difficult to meet the high accuracy requirements of certain projects, such as precision agriculture and crop yield estimation. Secondly, most studies do not strictly distinguish these two dates due to the short span from heading to flowering (usually approximately a week). In some studies [13,19,28], the VI_max method is used to extract the heading date, while in others [1,16,17], this method is used to extract the flowering date, which causes uncertainty in the extraction results of the heading date and the flowering date. Thirdly, the climate conditions (especially temperature) have a great impact on the growth and development of winter wheat. When the temperature is higher than 0 °C for consecutive days, winter wheat will begin to turn green and then develops to a specific growth stage when the accumulated temperature meets the requirements. However, the methods mentioned above are almost only based on remote sensing vegetation index data, so the advantages of meteorological data that is closely related to phenology have not been fully explored. Moreover, those methods can only realize monitoring and not forecasting.

Given that the growth of winter wheat is mainly driven by temperature of irrigated fields [28,29], and in a certain region, the accumulated temperature required for its development from green-up date to specific growth stages, such as heading and flowering, is relatively fixed, this paper aims to combine the satellite-derived green-up date and daily temperature data to develop a simple and accurate method to monitor and forecast the heading and flowering dates of winter wheat.

2. Study Area and Data

2.1. Study Area

The North China Plain (NCP) is located in northern and eastern China (110.4°E–122.7°E, 31.4°N–42.6°N) and lies in the lower reaches of the Yellow River. It is one of the three great plains in China, covering approximately 310,000 square kilometers (Figure 1). The NCP is an important grain production base in China. The main crops are winter wheat and summer maize. According to the National Bureau of Statistics of China [30], in 2018, the sown area of winter wheat in the NCP was 12.25 million ha, accounting for 53.9% of China, and the yield of winter wheat in the NCP was 75.76 million tons, accounting for 60.6% of China. The NCP is a typical area for planting winter wheat in China.

2.2. Data

2.2.1. Remote Sensing Data

The moderate-resolution imaging spectroradiometer (MODIS) MOD09Q1 product data from 2017 to 2019 were used to extract the planting area and the green-up date of winter wheat. The time resolution is 8 days, and the spatial resolution is 250 m. The widely used normalized difference vegetation index (NDVI), which is sensitive to vegetation signals in initial growth stage of crops [31], was calculated based on Equation (1):

$$\mathrm{NDVI}=\frac{{\rho }_{NIR}-{\rho }_R}{{\rho }_{NIR}+{\rho }_R}$$

(1)

where ${\rho }_{NIR}$ and ${\rho }_R$ represent the ground reflectance of the near-infrared band and red band, respectively.

Since noises existed in the NDVI time-series data because of the effects of clouds and atmospheric conditions, the change-weight filtering method [32], which can better preserve crop phenology characteristics, was used to denoise and smooth the NDVI time-series data. After that, the NDVI time-series data with good quality were prepared to extract the planting area and the green-up date of winter wheat.

2.2.2. Phenology Observation Data

The phenology observation data were acquired from the China Meteorological Administration (CMA). The dataset contains the field records of crop growth and development status since September 1991. The specific contents include crop variety, the name of the crop development stage and its date, the anomaly of the development stage, etc. In this paper, winter wheat observation stations recording both the green-up date and the heading date or the green-up date and the flowering date in NCP from 2017 to 2019 were selected. The stations with strong noise in the NDVI time-series were excluded to avoid interfering with the analysis. Finally, a total of 49 stations were used (Figure 1), and the number of available observation samples for each year is shown in

Table 1. Number of available observation samples.

Year	Number of Observation Samples
	Heading Date	Flowering Date
2017	37	36
2018	43	43
2019	48	48
Total	128	127

.

2.2.3. Meteorological Data

The meteorological data were derived from the dataset of CFSV2 (NCEP Climate Forecast System Version 2) on the Google Earth Engine (GEE) platform, which contains a variety of meteorological data worldwide since 1979, such as temperature, precipitation, and radiation [33]. The CFSV2 dataset is a common used dataset in related researches and its accuracy has been proven to be relatively high [34]. The time resolution of this dataset is 6 h, and the spatial resolution is 0.2°. Because the winter wheat in NCP is irrigated, the growth and development of winter wheat in NCP is mainly driven by air temperature, especially the daily minimum air temperature (denoted T_min) and the daily average air temperature (denoted T_mean). Therefore, the minimum temperature and the mean temperature of the CFSV2 dataset were used. Given the time resolution is 6 h, we took the minimum value of the four minimum temperatures within each day as the daily minimum temperature and the average value of the four mean temperatures within each day as the daily mean temperature.

2.2.4. Winter Wheat Map

The winter wheat map in NCP (Figure 1) was produced by a double threshold classification method based on the NDVI time-series data in 2019 (see details in Appendix A). The overall classification accuracy is 90.3%, producer accuracy is 86.8%, user accuracy is 96.6%, and Kappa coefficient is 0.80. Since we focus on the estimation of heading and flowering dates of winter wheat, it is important to ensure that the mapped pixels correspond to winter wheat though some winter wheat pixels may be not mapped, that is, we pursue higher user accuracy for the winter wheat map. In addition, the total planting area of winter wheat in the study area only changed less than 1% from 2017 to 2019 (http://www.stats.gov.cn), indicating that the distribution of winter wheat was quite stable during this period. Therefore, we used the same winter wheat map in 2019 for the three years of 2017, 2018, and 2019.

3. Methodology

3.1. The Proposed Accumulated Temperature Method

The accumulated temperature method (ATM) proposed in this study assumes that the growth and development of winter wheat is mainly driven by temperature (note that there are irrigation facilities in the study area, so the growth of winter wheat is not limited by water), and the accumulated value of the effective temperature required for its development from green-up date to specific growth stages (such as heading and flowering) is relatively fixed in a certain region. Therefore, the median accumulated temperature of observation samples is taken as the thermal requirement of the area. After that, the heading date and the flowering date of winter wheat can be estimated by combining the green-up date retrieved by remote sensing and daily temperature data.

3.1.1. Calculation of the Thermal Indices

The effective temperature (ET) and accumulated effective temperature (AET) for the heading/flowering date of winter wheat are calculated with Equations (2) and (3), respectively.

$${\mathrm{ET}}_i=\left\{\begin{array}{ll}{T}_i-T_{base},& T_i\ge T_{base}\\ 0,& T_i<T_{base}\end{array}\right.$$

(2)

$$\mathrm{AET}={\displaystyle \sum }_{i=g}^{k}E{T}_i$$

(3)

where $T_i$ represents the minimum or mean air temperature of the $i$th day, $T_{base}$ represents the baseline temperature, ${\mathrm{ET}}_i$represents the effective temperature of the $\mathrm{i}$th day, $g$ represents the green-up date, and $k$ represents the heading or flowering date.

Before calculating the ET and the AET, we need to determine the T_base. Considering the differences in winter wheat varieties and climate conditions of subregions in NCP, the T_base was calculated pixel by pixel. In this study, T_base is the average value of daily temperature during a given period before the green-up date. Two temperature indices, T_min and T_mean, and two periods, from October 1st last year (the sowing time of winter wheat) to the green-up date this year (denoted P₁) and one month before the green-up date (denoted P₂), were tested for examination. Based on the four combinations of two temperature indices and two periods for each phenological event (heading or flowering), we calculated the AET for each combination and named them as T_min_P₁, T_min_P₂, T_mean_P₁ and T_mean_P₂, respectively. We assume that a lower variation in the AETs among different observation samples indicates a stable T_base for calculation of the corresponding AET. Therefore, the coefficient of variation (CV) was used as a quantitative evaluation index to select the optimal T_base. Specifically, the T_base with a smaller CV in its calculated AETs was selected as the optimal one.

After the T_base is determined, the AET for the heading or flowering date can be calculated with the phenology observation samples. According to the definition of abnormal values in the boxplot, the abnormal samples whose AET was in excess of 1.5 times the interquartile range (IQR) were eliminated. Then, the median AET of the remaining samples was taken as the thermal requirement of the area. Moreover, to obtain a relatively stable thermal requirement, T_base was the average value of T_base for three years and the thermal requirement was also determined by samples of three years instead of year by year.

3.1.2. Determination of the Heading Date and Flowering Date

The heading date and flowering date can be determined by combining the daily air temperature data, the thermal requirement of the area, and the green-up date. The green-up date was extracted by the dynamic threshold method [12] based on NDVI time-series data, which was widely used in relevant studies [35,36,37]. According to the previous studies [38,39,40], the 20% threshold was used to extract the green-up date of winter wheat.

The heading/flowering date is the date whose accumulated value of ET from the retrieved green-up date exceeds the AET_median of heading/flowering for the first time. The calculation method is shown in Equation (4).

$$\mathrm{Pheno}=d+\mathrm{days}\left(AE{T}_{median}\right)$$

(4)

where $\mathrm{Pheno}$ is the heading/flowering date, $d$ represents the green-up date retrieved by remote sensing, $AE{T}_{median}$ represents the median $AET$ of effective observation samples in the area for the heading/flowering date, and $\mathrm{days}\left(\right)$ represents the number of days required for the accumulated value of $ET$ to exceed $AE{T}_{median}$ for the first time since the $d$th day.

3.2. Assessment of Monitoring Accuracy

Due to the large difference in the spatial scale between station observation and remote sensing, a given phenology observation station may be not directly corresponding to the pixel in which it is located. Therefore, we made a 3 km buffer for each phenology station. If there are winter wheat pixels within the buffer, the estimated phenology dates (green-up date, heading date, and flowering date) in those winter wheat pixels are averaged to represent the estimated phenology date for the phenology station.

Four indicators were used for accuracy assessment, including the coefficient of determination (R²), regression coefficient (a), BIAS, and root mean square error (RMSE). R² and a reflect the consistency between the estimated result and the observation data. The closer the R² and a are to 1, the higher the consistency between them. BIAS is defined as the number of days that the estimated result deviates from the observation data. If the deviation is greater than 0, the estimated result is later than the ground observation. If the deviation is less than 0, the estimated result is earlier than the ground observation. A smaller absolute value of BIAS indicates a higher accuracy of the estimated results. RMSE is the average error between the estimated result and the observation data. A smaller RMSE means a higher estimating accuracy. The calculations of the four indicators are shown in Equations (5)–(8), respectively.

$${\mathrm{R}}^2=\frac{cov{\left(\widehat{Y}, Y\right)}^2}{var\left(\widehat{Y}\right)var\left(Y\right)}$$

(5)

$$y=ax+b$$

(6)

$$\mathrm{BIAS}=\frac{{{\displaystyle \sum }}_{i=1}^N\left({\widehat{Y}}_i-Y_i\right)}{N}$$

(7)

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{\left({\widehat{Y}}_i-Y_i\right)}^2}{N}}$$

(8)

where ${\widehat{Y}}_i$ is the estimated result of the $i$th sample, $Y_i$ is the corresponding ground observation value, $N$ is the number of samples, $cov\left(\widehat{Y}, Y\right)$ represents the covariance between the estimated results and the ground observation values, and $var\left(\widehat{Y}\right)$ and $var\left(Y\right)$ represent the variance of the estimated results and the ground observation values, respectively. $y=ax+b$ is the linear regression model of the estimated results and the ground observation values in which the estimated result is the dependent variable and the ground observation value is the independent variable.

In addition, the VI_max method [13], which is usually used to extract the heading–flowering date of winter wheat in current researches, was selected as the comparison method in this study. The VI_max method considers that the heading–flowering of crops is the transformation phase from vegetative growth to reproductive growth, and leaves begin to wither and die after the heading–flowering phase. Therefore, the corresponding date of the maximum vegetation index in the growing season of crops is assumed as the heading–flowering date [13]. It should be noted that the VI_max method cannot distinguish between the heading date and the flowering date. The extracted date is roughly in the period from heading to flowering, so the estimated result is compared with both the observed heading date and the observed flowering date (Section 4.2.2).

3.3. Assessment of Forecasting Accuracy

To test whether the proposed method has good ability to forecast the heading and flowering dates of winter wheat, we used the thermal requirement derived from previous years to forecast the phenology in subsequent years and evaluated the accuracy of the forecasting results. The experiment was divided into two parts. The first was to forecast the heading and flowering dates in 2018 and 2019 using the thermal requirement derived from the observation samples in 2017. The second was to forecast the heading and flowering dates in 2019 based on observation samples in both 2017 and 2018. The spatial scale match method and accuracy indicators are the same as the assessment of monitoring accuracy.

4. Results

4.1. Performance of Different T_base

For the heading date, the CV based on different T_base ranges from 0.14 to 0.41 (Figure 2a), and for the flowering date, the CV based on different T_base ranges from 0.13 to 0.35 (Figure 2b). For both the heading date and the flowering date, the temperature indices based on P₂ outperform the temperature indices based on P₁ since the AETs based on P₂ are much more concentrated in the frequency distribution and their CVs are much smaller. Similarly, comparing T_mean with T_min, the AETs based on T_mean are more concentrated and their CVs are much smaller. Therefore, T_mean_P₂ (the average value of daily mean temperature during the month before the green-up date) is selected as the optimal T_base, and the corresponding AET is set as 527.4 °C for the heading date and 628.7 °C for the flowering date (Figure 2). The optimal T_base in NCP is from about 0 to 6 °C (Figure 3a). It gradually increases from the south to the north, which is similar to the spatial distribution pattern of green-up date of winter wheat in NCP (Figure 3b).

4.2. Evaluation of Monitoring Accuracy

4.2.1. The Spatial Distribution of Estimated Heading and Flowering Dates

The spatial distribution patterns of both the heading and flowering dates estimated by the ATM method among different years are similar (Figure 4a–f), while the spatial distribution patterns of the estimated heading–flowering date by the VI_max method vary largely over time (Figure 4g–i). In addition, the heading and flowering dates estimated by the ATM method show a higher spatial continuity (Figure 4(b₁,e₁)) and a smoother transition from the south to the north (Figure 4a–f) compared with the VI_max method. The estimated heading–flowering date with the VI_max method show a large variation even in a small region (Figure 4(h₁)) and their spatial transitions from the south to the north are abrupt (Figure 4g–i).

4.2.2. Monitoring Accuracy of Estimated Phenology

For each year, the monitoring accuracy for the heading/flowering date based on the ATM method is much higher than that based on the VI_max method (Figure 5 and Figure 6). Compared with the VI_max method, the R² for the estimated heading (flowering) date by the ATM method increases by 0.25 to 0.34 (0.21 to 0.23) among different years, the RMSE decreases by 2.09 to 6.97 d (1.61 to 8.08 d), the a is closer to 1, and the BIAS is closer to 0.

4.3. Evaluation of Forecasting Accuracy

When only using samples in 2017, the R² for forecasted phenology ranges from 0.54 to 0.69, a ranges from 0.79 to 1.18, BIAS ranges from −2.02 to −2.94 d, and RMSE ranges from 6.13 to 6.38 d (

Table 2. Accuracy of the forecasted phenology by the ATM method.

Phenology	Year	R2	a	BIAS/d	RMSE/d
	Based on samples in 2017
Heading date	2018	0.69	1.18	−2.79	6.38
	2019	0.62	0.79	−2.94	6.22
Flowering date	2018	0.54	1.06	−2.42	6.36
	2019	0.57	0.92	−2.02	6.13
	Based on samples in 2017 and 2018
Heading date	2019	0.60	0.77	−0.65	5.62
Flowering date	2019	0.58	0.93	−0.21	5.78

). The forecasting accuracy based on samples of the previous year is lower than that based on samples of the current year (i.e., monitoring accuracy), but it is still at a relatively high level. Compared with the thermal requirement derived from samples of the previous single year, the forecasting accuracy based on the thermal requirement derived from samples of the previous two years is generally higher. When using samples in 2017 and 2018 together, the R² and a have no obvious change, but the RMSE decreases by 0.60 d and 0.35 d for the heading date and flowering date, respectively, and the BIAS is much closer to 0.

5. Discussion

5.1. Advantages of the Proposed Method

Compared with the widely used VI_max method, the estimated phenology with the ATM method is more accurate. The visual qualitative analysis shows that the spatial distribution patterns of the estimated heading and flowering dates with the ATM method are similar among different years and have a higher spatial continuity. By contrast, the spatial inconsistency in the estimated heading–flowering date with the VI_max method indicates that it is more susceptible to the noise in the VI time-series data. Normally, the phenology will not change greatly during the adjacent years [41,42]. It is also relatively stable in a small region and is consistent with the climate gradient [43,44]. Therefore, the estimated phenology with the ATM method is more consistent with the reality. Furthermore, the quantitative analysis shows that the RMSE of the estimated heading and flowering dates with the ATM method are approximately 5–6 d, and the BIAS values are within ±1 d. By contrast, there is a relatively large difference between the heading–flowering dates estimated with the VI_max method and the ground observations. The RMSE varies from 8 to 13 d, and the BIAS is between −7 and 8 d. It should be mentioned that the VI_max method in its original reference [13] was based on the enhanced vegetation index (EVI) time-series data. Due to the anti-saturation ability of EVI, for the VI_max method, it may be better to use EVI to extract the heading–flowering date. However, it will not make a big difference between the results based on NDVI and those based on EVI. The RMSE was approximately 10 d even if the EVI was used [13], which was close to our experiment (Figure 5 and Figure 6).

Considering the differences between remote sensing and in situ meteorological data, we compared the RMSE based on AET (in situ), AET (remote sensing), and VI_max (see details in Appendix B). The results show that the RMSE based on AET (in situ) was very close to that based on AET (remote sensing) but a little lower. This may be because the resolution of the remote sensing meteorological dataset is relatively coarse (0.2°), and the accuracy is lower than that of the ground observation stations. It also indicated both the remote sensing and in situ meteorological data can be effectively used in the ATM method and generate results with a high accuracy.

The ATM method can well distinguish between the physiological heading date and flowering date and has an obvious advantage in distinguishing phenological events which are close in time. By contrast, the VI_max method can only extract the date corresponding to the maximum vegetation index in the growing season, which is essentially not a clear crop phenological event and does not have physiological significance. It can hardly distinguish the phenological events (such as the heading and the flowering), which are close in time.

Additionally, the ATM method can realize not only monitoring but also forecasting. Compared with the monitoring accuracy, the forecasting accuracy decreases slightly but is still at a high level (RMSE is approximately 6 d). When the green-up date is known, meteorological observations or forecast data can be used to estimate the phenology. Therefore, combination with the short-term meteorological forecast data can realize forecasting of the heading and flowering dates of winter wheat, while the existing methods based on remote sensing data can only realize monitoring.

5.2. Limits and Future Improvements of the Proposed Method

Although the ATM method can produce more accurate results than the VI_max method, it relies on extra meteorological data in addition to remote sensing vegetation index data. Moreover, the spatial resolution of meteorological data is relatively coarse compared with remote sensing data. However, with the launch of meteorological satellites in recent years, daily land surface temperature (LST) products with finer spatial resolutions have become increasingly abundant (e.g., MOD11A1, FY3B, etc.). Considering there are still differences between the surface temperature and air temperature as well as the surface temperature products and ground observed temperatures [45], whether these products could be used to monitor and forecast the heading date and the flowering date of winter wheat is needed to be further studied.

In addition, when using the ATM method to monitor and forecast crop phenology, daily air temperature data and the thermal requirements for crop growth and development are needed in advance according to the historical phenology observation data or phenology calendar in the region. Compared with the VI_max method, the ATM method requires some prior knowledge. In this study, the winter wheat in NCP is irrigated, while the accuracy of this method when it is applied to the rain-fed croplands has not been investigated. Moreover, whether this method can still be effective in years when the hydrothermal conditions are greatly different from the average level (e.g., sudden natural disasters such as drought, cold damage, and floods) needs to be further verified.

6. Conclusions

In this study, an accumulated temperature method (ATM) for monitoring and forecasting the heading and flowering dates of winter wheat was proposed from the perspective of the thermal requirements for the growth of crops. The method combined the green-up date derived from remote sensing data and daily temperature data to estimate and forecast the phenology. First, the green-up date extracted by the dynamic threshold method was taken as the starting point of temperature accumulation. Then, the accumulated effective temperature from the green-up date to the heading or flowering date for each phenology observation sample was calculated. After that, the median value was selected as the thermal requirement in the study area. Finally, the daily temperature data was used to estimate the heading and flowering dates of winter wheat. The phenology estimated by the ATM method is relatively stable among adjacent years, and the spatial distribution shows good continuity and smooth transition among different local areas. The monitoring accuracy of the ATM method is high, and it can effectively distinguish between the heading date and the flowering date which are close in time. In addition, the method also performs well when used to forecast the heading and flowering dates of winter wheat with the short-term meteorological forecast data.