Determination of Key Phenological Phases of Winter Wheat Based on the Time-Weighted Dynamic Time Warping Algorithm and MODIS Time-Series Data

: Accurate determination of phenological information of crops is essential for ﬁeld management and decision-making. Remote sensing time-series data are widely used for extracting phenological phases. Existing methods mainly extract phenological phases directly from individual remote sensing time-series, which are easily affected by clouds, noise, and mixed pixels. This paper proposes a novel method of phenological phase extraction based on the time-weighted dynamic time warping (TWDTW) algorithm using MODIS Normalized Difference Vegetation Index (NDVI) 5-day time-series data with a spatial resolution of 500 m. Firstly, based on the phenological differences between winter wheat and other land cover types, winter wheat distribution is extracted using the TWDTW classiﬁcation method, and the results show that the overall classiﬁcation accuracy and Kappa coefﬁcient reach 94.74% and 0.90, respectively. Then, we extract the pure winter-wheat pixels using a method based on the coefﬁcient of variation, and use these pixels to generate the average phenological curve. Next, the difference between each winter-wheat phenological curve and the average winter-wheat phenological curve is quantitatively calculated using the TWDTW algorithm. Finally, the key phenological phases of winter wheat in the study area, namely, the green-up date (GUD), heading date (HD), and maturity date (MD), are determined. The results show that the phenological phase extraction using the TWDTW algorithm has high accuracy. By veriﬁcation using phenological station data from the Meteorological Data Sharing Service System of China, the root mean square errors (RMSEs) of the GUD, HD, and MD are found to be 9.76, 5.72, and 6.98 days, respectively. Additionally, the method proposed in this article is shown to have a better extraction performance compared with several other methods. Furthermore, it is shown that, in Hebei Province, the GUD, HD, and MD are mainly affected by latitude and accumulated temperature. As the latitude increases from south to north, the GUD, HD, and MD are delayed, and for each 1 ◦ increment in latitude, the GUD, HD, and MD are delayed by 4.84, 5.79, and 6.61 days, respectively. The higher the accumulated temperature, the earlier the phenological phases occur. However, latitude and accumulated temperature have little effect on the length of the phenological phases. Additionally, the lengths of time between GUD and HD, HD and MD, and GUD and MD are stable at 46, 41, and 87 days, respectively. Overall, the proposed TWDTW method can accurately determine the key phenological phases of winter wheat at a regional scale using remote sensing time-series data.


Introduction
Wheat is one of the most important food crops in the global market. Winter wheat accounts for about 80% of global wheat production [1]. Crop phenology is an important link between crop ecosystems and yield estimation. Accurate acquisition of winter-wheat phenological information is of great significance for field management, decision-making, and yield prediction. Traditionally, winter-wheat phenological information is extracted using costly and time-consuming methods, such as manual recording and meteorological station observations, which are not suitable for investigating large areas [2]. It is difficult to accurately extract the phenological phases of winter wheat over a large area, due to the influence of latitude and climatic factors.
Remote sensing platforms can obtain measurements quickly and over large areas, and remote sensing data are widely used in various kinds of studies. For example, such data are convenient for obtaining long-term vegetation growth information. Long-term observational data of vegetation can help determine the response of plants to climate change, and such data are widely used in vegetation phenology research-for example, time-series of the Normalized Difference Vegetation Index (NDVI) [3][4][5][6][7] and Enhanced Vegetation Index (EVI) [8,9]. Existing methods for extracting phenological information from long-term remote sensing time-series mainly include threshold methods and change detection methods [2,10,11].
Threshold methods are a simple and easy-to-implement way to extract phenological indicators from remote sensing time-series. The earliest threshold method is the fixed threshold method, which establishes a fixed threshold according to the phenological characteristics of vegetation. However, this method has a shortcoming, namely, that the threshold is difficult to set, due to the inter-annual variability of vegetation and climate change. To address this limitation, the dynamic threshold method was developed, which is related to the seasonal variation of vegetation. Vrieling et al. [12] used a hyperbolic tangent model to better fit the upper envelope of observations and define the start of season (SOS) as the moment when the NDVI reaches 20% of its maximum annual amplitude and define the peak season (PS) as the time when the NDVI reaches 90% of its maximum annual amplitude. Additionally, Gan et al. [13] used six methods, including the relative threshold at 10%, 20%, or 50% of the amplitude of the vegetation index (RT10, RT20, and RT50, respectively), to detect the green-up date (GUD) of winter wheat for the period of 2009-2013 in the Huanghuaihai winter-wheat region of China. Compared with the fixed threshold method, the dynamic threshold method fully considers the characteristics of the vegetation index curve of the target crop and sets the threshold according to the specific conditions of the research area, which can eliminate the influence of soil background and other crops. However, when the crop phenology is extracted for complex terrain, mixed pixel phenomena will affect the accuracy of crop phenology monitoring. The change detection method involves fitting the vegetation index curve and then determining the change point of the curve in combination with the vegetation phenological characteristics and the change of the fitting curve. Before phenological information is extracted, long remote sensing time-series data are smoothed. Data smoothing can eliminate noise and reconstruct a more representative vegetation curve to provide key parameters for phenological monitoring. There have been numerous studies on vegetation curve smoothing methods, and the most widely used methods include logistic fitting [14,15], and its improved algorithm [16][17][18], the polynomial fitting algorithm [19][20][21][22][23], and the asymmetric Gaussian algorithm [24]. Typical change detection methods include the curvature method [14], and the slope method [25]. For example, Zhang et al. [14] first used the segmental logistic function to fit the Moderate Resolution Imaging Spectroradiometer (MODIS) vegetation index and then used the phenological period of vegetation corresponding to the phenological curvature to calculate the extreme point of curvature change as the GUD, maturity date (MD), decay date, and dormancy date of vegetation growth. Furthermore, Lu et al. [25] determined the SOS of winter wheat in the North China Plain as the time corresponding to the maximum derivative of the NDVI time-series curve. All the aforementioned methods use the time-series of each pixel to extract the phenology and are influenced by clouds, noise, and mixed pixels. Additionally, existing research extracts and analyzes only one phenological phase and lacks a collaborative analysis of multiple phenological phases.
The time-weighted dynamic time warping (TWDTW) algorithm is an improvement of the Dynamic Time Warping (DTW) algorithm [26]. The TWDTW algorithm was recently applied to remotely sensed time-series data [27][28][29][30][31]. TWDTW algorithm works by comparing the similarity between two time-series and finds their optimal ailgnment, resulting in dissimilarity measure. Studies have shown that the TWDTW algorithm is more robust in crop and forest classification than the commonly used machine learning algorithms (e.g., support vector machine (SVM), random forest (RF)), and can still achieve higher classification accuracy with limited training samples [32][33][34]. However, the TWDTW algorithm has never been used for crop phenology extraction. In the extraction of phenology, for the same crop, due to the difference in geographical locations, climatic factors, and field management, the shape of the crop phenological curve will change. The TWDTW algorithm can provide a matching relationship when calculating the similarity of two curves to obtain the time alignment relationship, and this relationship can be used to calculate the difference in phenology. Thus, the accurate phenological phase can be determined. Overall, TWDTW is used both for classification and phenology extraction in this study.
In this study, we propose a novel method for phenological phase extraction based on the TWDTW algorithm. We firstly extracted winter wheat distribution using the TWDTW classification method [32][33][34]. Then, we extracted the pure winter-wheat pixels using a method based on the coefficient of variation and then used these pixels to generate the average phenological curve. Furthermore, the GUD, heading date (HD), and MD were defined based on the average phenological curve. Next, waveform processing was performed on each NDVI time-series curve. Finally, the difference between each NDVI phenological curve and the average phenological curve was extracted using the TWDTW algorithm to obtain the phenological phases.

Study Area
The study area is located in Hebei Province in Northern China, between 36 • 05 and 42 • 40 N and between 113 • 27 and 119 • 50 E ( Figure 1). Hebei Province contains plateaus, mountains, hills, and plains, with the latter accounting for 41.2% of the total area of the province. Hebei Province has a temperate continental monsoon climate, with cold and dry winters, dry and windy springs, hot and rainy summers, and clear and cool autumns. The province's average annual rainfall is 350-815 mm, and its annual average temperature is (-0.3)-14.0 • C. The main local cultivation system is a double-season rotation, with winter wheat-corn rotation being the most common; other crops include cotton, spring corn, millet, and peanuts.

Remote Sensing Data and Preprocessing
Remote sensing data were obtained from the Terra satellite daily surface reflectivity product MOD09GA, which is developed by the National Aeronautics and Space Administration (NASA) MODIS land product group (https://ladsweb.modaps.eosdis.nasa.gov, accessed on 26 June 2020). The product includes daily surface reflectance data from channels 1-7 with a spatial resolution of 500 m. The MODIS NDVI daily product was obtained from the red and near-infrared channels, and the maximum value of each five-day period was then taken to obtain the MODIS NDVI five-day composite products. A total of 31 NDVI composites from 01 February to 30 June 2018 were used in this study, which covers the key phenological growth stages of winter wheat, such as GUD, HD, and MD. To reduce the influence of noise on the NDVI time-series data, a Savitzky-Golay (S-G) [13,35,36] filter was applied to the whole time-series. In this study, the MODIS data were mainly used to determine the winter-wheat distribution and phenological phases.

Remote Sensing Data and Preprocessing
Remote sensing data were obtained from the Terra satellite daily surface reflectivity product MOD09GA, which is developed by the National Aeronautics and Space Administration (NASA) MODIS land product group (https://ladsweb.modaps.eosdis.nasa.gov, accessed on 26 June 2020). The product includes daily surface reflectance data from channels 1-7 with a spatial resolution of 500 m. The MODIS NDVI daily product was obtained from the red and near-infrared channels, and the maximum value of each five-day period was then taken to obtain the MODIS NDVI five-day composite products. A total of 31 NDVI composites from 1 February to 30 June 2018 were used in this study, which covers the key phenological growth stages of winter wheat, such as GUD, HD, and MD. To reduce the influence of noise on the NDVI time-series data, a Savitzky-Golay (S-G) [13,35,36] filter was applied to the whole time-series. In this study, the MODIS data were mainly used to determine the winter-wheat distribution and phenological phases.

Phenological Monitoring Station Data and Field Data
The phenological monitoring station data were obtained from the Meteorological Data Sharing Service System of China (http://data.cma.cn, accessed on 4 August 2020). This data set contains field records of crop growth and development since September 1991. The data include crop varieties, names, and dates of crop development stages, abnormal development stages, etc. Data from 15 phenological monitoring stations in Hebei Province were used (Figure 1b), each of which recorded the GUD, HD, and MD. The ground GUD was recorded by observers when more than half of the winter-wheat leaves had begun to turn green and reached a length of 1-2 cm in a wheat field near the station; the ground HD was recorded by observers when more than half of the top of the ear (excluding the awns) was exposed from the leaf sheath in a wheat field near the station; and the ground MD was recorded by observers when the endosperm of more than 50% of grains was waxy, and the ear and the node under the ear had turned yellow in a wheat field near the station. These data were used to verify the results for the phenological phases extracted

Phenological Monitoring Station Data and Field Data
The phenological monitoring station data were obtained from the Meteorological Data Sharing Service System of China (http://data.cma.cn, accessed on 4 August 2020). This data set contains field records of crop growth and development since September 1991. The data include crop varieties, names, and dates of crop development stages, abnormal development stages, etc. Data from 15 phenological monitoring stations in Hebei Province were used (Figure 1b), each of which recorded the GUD, HD, and MD. The ground GUD was recorded by observers when more than half of the winter-wheat leaves had begun to turn green and reached a length of 1-2 cm in a wheat field near the station; the ground HD was recorded by observers when more than half of the top of the ear (excluding the awns) was exposed from the leaf sheath in a wheat field near the station; and the ground MD was recorded by observers when the endosperm of more than 50% of grains was waxy, and the ear and the node under the ear had turned yellow in a wheat field near the station. These data were used to verify the results for the phenological phases extracted using the TWDTW method. Additionally, we collected ground samples from field surveys in 2018. The location of each winter wheat sample was recorded using a mobile GPS device. Then, we added some non-winter-wheat samples (e.g., water, forest grass, town) through visual interpretation on the high-resolution image of Google Earth. Finally, 263 winter-wheat samples and 223 non-winter-wheat samples were collected. The samples were divided into two parts-one part used for training and classification, another part used for accuracy evaluation. (Figure 1b).

Temperature Data
The temperature data were based on European Space Agency (ESA) EAR5 data (https://cds.climate.copernicus.eu/, accessed on 16 August 2020). Data were acquired for latitudes between 35 • and 45 • N and longitudes between 112 • and 120 • E, covering the whole of Hebei Province. The grid intervals of longitude and latitude were both 0.1 • . The EAR5 data provide one temperature measurement per hour, that is, 24 data sets every day. Data were obtained from 01 February to 30 June 2018, that is, for a total of 150 days. In Remote Sens. 2021, 13, 1836 5 of 18 this study, the temperature data were mainly used to analyze the effects of accumulated temperature on the winter-wheat phenological phases.

Extraction of Winter Wheat Distribution by TWDTW Classification Method
The land cover sample points and MODIS NDVI time-series data were used to construct the NDVI time-series curve. The waveform of the time-series curve reflects the growth and changes of land cover types. According to the phenological characteristics of different land cover types, the winter wheat distribution was extracted using the TWDTW classification method.

Phenological Characteristics of Land Cover
Different land cover types have different phenological patterns. Therefore, it is very important to understand the changes in the phenological characteristics of different land cover types [37]. During the growth period of winter wheat, the NDVI time-series curves of several typical land cover types are different from those of winter wheat. As shown in Figure 2, the NDVI of winter wheat is greater than the NDVI of non-vegetation land cover types (e.g., water, town). The NDVI of winter wheat gradually increased from March, peaked in May, then gradually decreased and reached the lowest value in mid-late June. In contrast, the NDVI of forest grass gradually increased from April and showed different phenological characteristics from winter wheat.
The temperature data were based on European Space Agency (ESA) EAR5 data (https://cds.climate.copernicus.eu/, accessed on 16 August 2020). Data were acquired for latitudes between 35° and 45° N and longitudes between 112° and 120° E, covering the whole of Hebei Province. The grid intervals of longitude and latitude were both 0.1°. The EAR5 data provide one temperature measurement per hour, that is, 24 data sets every day. Data were obtained from 1 February to 30 June 2018, that is, for a total of 150 days. In this study, the temperature data were mainly used to analyze the effects of accumulated temperature on the winter-wheat phenological phases.

Extraction of Winter Wheat Distribution by TWDTW Classification Method
The land cover sample points and MODIS NDVI time-series data were used to construct the NDVI time-series curve. The waveform of the time-series curve reflects the growth and changes of land cover types. According to the phenological characteristics of different land cover types, the winter wheat distribution was extracted using the TWDTW classification method.

Phenological Characteristics of Land Cover
Different land cover types have different phenological patterns. Therefore, it is very important to understand the changes in the phenological characteristics of different land cover types [37]. During the growth period of winter wheat, the NDVI time-series curves of several typical land cover types are different from those of winter wheat. As shown in Figure 2, the NDVI of winter wheat is greater than the NDVI of non-vegetation land cover types (e.g., water, town). The NDVI of winter wheat gradually increased from March, peaked in May, then gradually decreased and reached the lowest value in mid-late June. In contrast, the NDVI of forest grass gradually increased from April and showed different phenological characteristics from winter wheat.

TWDTW Classification Method
The TWDTW algorithm includes a time weighting factor based on dynamic warping to find the path with the minimum cumulative distance. When calculating the degree of similarity between two sequences, it is important to consider not only the numerical size of the matching points, but also the period between the matching points to avoid a severe mismatch between sequences [32]. Suppose that there are two time-series curves: X = (X1, X2, X3, . . . , Xm), and Y = (Y1, Y2, Y3, . . . , Yn). First, the distance matrix of the two time-series D m×n [Equation (1)] ( Figure 3b) is constructed, where each element in the matrix is defined as D i,j = w i,j * X i − Y j , where w i,j is the time weighting factor following Remote Sens. 2021, 13, 1836 6 of 18 a logistic function, as given in Equation (2). Secondly, the cumulative distance matrix C m×n [Equation (3)] (Figure 3c) is calculated.
where α and β represent the steepness and median value of the time-series [32]. Thirdly, the minimum distance between the two curves and the corresponding matching relationship is determined using the cumulative distance matrix C m×n (Figure 3c). According to this matching relationship, the time alignment relationship of the two curves is obtained (Figure 3d).
to find the path with the minimum cumulative distance. When calculating the degree of similarity between two sequences, it is important to consider not only the numerical size of the matching points, but also the period between the matching points to avoid a severe mismatch between sequences [32]. Suppose that there are two time-series curves: X = (X1, X2, X3,…, Xm), and Y = (Y1, Y2, Y3,…, Yn). First, the distance matrix of the two time-series Figure 3b) is constructed, where each element in the matrix is defined as , = , * − , where , is the time weighting factor following a logistic function, as given in Equation (2). Secondly, the cumulative distance matrix where α and β represent the steepness and median value of the time-series [32]. Thirdly, the minimum distance between the two curves and the corresponding matching relationship is determined using the cumulative distance matrix × (Figure 3c). According to this matching relationship, the time alignment relationship of the two curves is obtained (Figure 3d).  Specifically, winter wheat in Hebei province was mapped using the following steps: (1) One hundred and forty-two winter-wheat samples and one hundred and sixteen nonwinter-wheat samples were randomly selected for training. One hundred and twentyone winter-wheat samples and one hundred and seven non-winter-wheat samples were selected for validation; (2) the average winter-wheat phenological curve was generated using the NDVI time-series curve of 142 winter-wheat samples; (3) the TWDTW distance between the average winter-wheat phenological curve and the phenological curve of the training samples were calculated, and the best threshold 1.5 is obtained; (4) the TWDTW distance between the average winter-wheat phenological curve and all the NDVI time-Remote Sens. 2021, 13, 1836 7 of 18 series curves was calculated, and the samples whose distance was less than the threshold were considered to be winter wheat.

Determination of Winter-Wheat Key Phenological Phases by TWDTW Algorithm
Based on the MODIS NDVI time-series data and the winter-wheat distribution, the key phenological phases of winter wheat were determined, including the following steps: (a) The selection of pure winter-wheat pixels; (b) the definition of the average GUD, HD, and MD; (c) the waveform adjustment; (d) the calculation of the difference between each NDVI phenological curve and the average phenological curve using the TWDTW algorithm; and (e) the determination of the key winter-wheat phenological phases and accuracy assessment (Figure 4).
winter-wheat samples were randomly selected for training. One hundred and twenty-one winter-wheat samples and one hundred and seven non-winter-wheat samples were selected for validation; (2) the average winter-wheat phenological curve was generated using the NDVI time-series curve of 142 winter-wheat samples; (3) the TWDTW distance between the average winter-wheat phenological curve and the phenological curve of the training samples were calculated, and the best threshold 1.5 is obtained; (4) the TWDTW distance between the average winter-wheat phenological curve and all the NDVI timeseries curves was calculated, and the samples whose distance was less than the threshold were considered to be winter wheat.

Determination of Winter-Wheat Key Phenological Phases by TWDTW Algorithm
Based on the MODIS NDVI time-series data and the winter-wheat distribution, the key phenological phases of winter wheat were determined, including the following steps: (a) The selection of pure winter-wheat pixels; (b) the definition of the average GUD, HD, and MD; (c) the waveform adjustment; (d) the calculation of the difference between each NDVI phenological curve and the average phenological curve using the TWDTW algorithm; and (e) the determination of the key winter-wheat phenological phases and accuracy assessment (Figure 4).

Selection of Pure Winter-Wheat Pixels
The accuracy of the average phenological curve directly affects the accuracy of the extraction of GUD, HD, and MD. In this study, large homogeneous areas of winter wheat were chosen to avoid mixed pixels, due to the coarser resolution of MODIS [38,39]. Using pure winter-wheat pixels to generate an average phenological curve can improve the accuracy of the extraction of GUD, HD, and MD. In this study, the coefficient of variation was used to obtain pure winter-wheat pixels [13]. The coefficient of variation, also known Remote Sens. 2021, 13, 1836 8 of 18 as the coefficient of dispersion, is defined as the ratio of the standard deviation to the mean. The steps in the process for the selection of pure winter-wheat pixels are as follows: (1) Establish a series of 3 × 3 windows with each winter-wheat pixel as the center. (2) Calculate the coefficient of variation of each of the nine pixels in each 3 × 3 window. As the time range of the MODIS NDVI time-series data used in this paper is from 1 February to 30 June 2018, and the time interval is five days, there are 31 images in total, and 31 coefficients of variation can be obtained for each pixel. (3) If the 31 coefficients of variation of a given pixel are all less than the threshold of 0.5%, the winter-wheat pixel is determined to be a pixel of pure winter-wheat.
Here, we argue that it is appropriate to set the size of the window to 3 × 3, considering that the winter-wheat planting area in Hebei Province is mainly in the plain area, and the wheat is planted continuously. If the window is too large, it is likely to be affected by mixed pixels, and no pure winter-wheat pixels can be found. Additionally, the threshold value of 0.5% was chosen-since, if the threshold value is too large, the influence of mixed pixels cannot be effectively removed; in contrast, if the threshold is too small, too few pixels of pure winter-wheat are extracted, which is not representative.

Definition of the Average GUD, HD, and MD
Based on the extracted pure winter-wheat pixels, the average winter-wheat phenological curve was obtained using the arithmetic average, which is used to define the average of GUD, HD, and MD. Since the average winter-wheat phenological curve is a time-series with a 5-day interval, a 6 th -order polynomial function was used to fit the 5-day interval NDVI time-series and reconstruct daily NDVI time-series; this approach has previously been used to extract crop phenology [19,[21][22][23]. Here, the GUD and MD were defined using the relative threshold method [13,24,35], and the HD was defined via the maximum method [35]. The definition of the average GUD, HD, and MD is as shown in Figure 5.
The accuracy of the average phenological curve directly affects the accuracy of the extraction of GUD, HD, and MD. In this study, large homogeneous areas of winter wheat were chosen to avoid mixed pixels, due to the coarser resolution of MODIS [38,39]. Using pure winter-wheat pixels to generate an average phenological curve can improve the accuracy of the extraction of GUD, HD, and MD. In this study, the coefficient of variation was used to obtain pure winter-wheat pixels [13]. The coefficient of variation, also known as the coefficient of dispersion, is defined as the ratio of the standard deviation to the mean. The steps in the process for the selection of pure winter-wheat pixels are as follows: (1) Establish a series of 3 × 3 windows with each winter-wheat pixel as the center. (2) Calculate the coefficient of variation of each of the nine pixels in each 3 × 3 window. As the time range of the MODIS NDVI time-series data used in this paper is from 01 February to 30 June 2018, and the time interval is five days, there are 31 images in total, and 31 coefficients of variation can be obtained for each pixel. (3) If the 31 coefficients of variation of a given pixel are all less than the threshold of 0.5%, the winter-wheat pixel is determined to be a pixel of pure winter-wheat.
Here, we argue that it is appropriate to set the size of the window to 3 × 3, considering that the winter-wheat planting area in Hebei Province is mainly in the plain area, and the wheat is planted continuously. If the window is too large, it is likely to be affected by mixed pixels, and no pure winter-wheat pixels can be found. Additionally, the threshold value of 0.5% was chosen-since, if the threshold value is too large, the influence of mixed pixels cannot be effectively removed; in contrast, if the threshold is too small, too few pixels of pure winter-wheat are extracted, which is not representative.

Definition of the Average GUD, HD, and MD
Based on the extracted pure winter-wheat pixels, the average winter-wheat phenological curve was obtained using the arithmetic average, which is used to define the average of GUD, HD, and MD. Since the average winter-wheat phenological curve is a time-series with a 5-day interval, a 6th-order polynomial function was used to fit the 5-day interval NDVI time-series and reconstruct daily NDVI time-series; this approach has previously been used to extract crop phenology [19,[21][22][23]. Here, the GUD and MD were defined using the relative threshold method [13,24,35], and the HD was defined via the maximum method [35]. The definition of the average GUD, HD, and MD is as shown in Figure 5. The average winter-wheat phenological curve is divided into two parts according to the maximum value of the curve, namely, SA and SB. The GUD was defined as the day when the NDVI first reached 10% of Amp1, and was obtained by subtracting NVDImin1 from NDVImax. The HD was defined as the day when the NDVI reached the maximum. The MD was defined as the day after the HD when the NDVI dropped to 10% of Amp2, and was obtained by subtracting NVDImin2 from NDVImax. Where NDVImax represents the maximum NDVI of the average winterwheat phenological curve, NDVImin1 and NDVImin2 represent the minimum NDVI of SA and SB, respectively. Amp1 and Amp2 represent the amplitude of SA and SB, respectively. The average winter-wheat phenological curve is divided into two parts according to the maximum value of the curve, namely, S A and S B . The GUD was defined as the day when the NDVI first reached 10% of Amp 1 , and was obtained by subtracting NVDI min1 from NDVI max . The HD was defined as the day when the NDVI reached the maximum. The MD was defined as the day after the HD when the NDVI dropped to 10% of Amp 2 , and was obtained by subtracting NVDI min2 from NDVI max . Where NDVI max represents the maximum NDVI of the average winter-wheat phenological curve, NDVI min1 and NDVI min2 represent the minimum NDVI of S A and S B , respectively. Amp1 and Amp2 represent the amplitude of S A and S B , respectively.

Waveform Adjustment
Before extracting the differences between each NDVI phenological curve and the average phenological curve, it is necessary to preprocess the phenological curve of winter wheat, including waveform adjustment and polynomial fitting. Due to factors, such as temperature, latitude, and longitude, and the NDVI value of winter wheat, may be different in different regions even for the same phenological phase. For example, in some regions, when winter wheat enters the HD, the NDVI reaches 0.8, while in other regions, the NDVI only reaches 0.7. Therefore, it was necessary to adjust the NDVI time-series curve to reduce the influence of different amplitudes on the two phenological curves. Details of the waveform adjustment method are given in Appendix A.

Calculation of the Difference between NDVI Phenological Curves and the Average Phenological Curve
According to the time alignment relationship by the TWDTW algorithm in Section 3.1.2, when the phenological phase of the reference curve is known, the phenological phase of the target curve can be obtained. Figure 6 shows an example of how the TWDTW algorithm was used to determine the phenological phase. Two NDVI phenological curves were used, namely, a reference curve and a target curve, where the key phenological phases (GUD, HD, and MD) of the reference curve (termed GUD R , HD R , and MD R , respectively) are known. Three cases are analyzed here. In the first case, one point on the reference curve matches one point on the target curve; HD R and the location of HD on the target curve (HD T ) belong to this situation; that is, there is only one line segment connecting HD R and HD T . This shows that, from the perspective of waveform similarity, HD T has similar characteristics to HD R , so HD T can be regarded as the HD of the target curve. In the second case, multiple points on the reference curve match one point on the target curve; GUD R and the location of GUD on the target curve (GUD T ) belong to this situation. GUD R and the three points on the left match the GUD T point. From the perspective of alignment, although GUD T is aligned with four points, GUD R only matches GUD T , so GUD T can be regarded as the GUD of the target curve. In the third case, one point on the reference matches multiple points on the target curve; this is the case for MD R and the location of MD on the target curve (MD T ). MD R is aligned with three points on the target curve. In this case, the average point is chosen as the MD of the target curve. After the phenological phase of the target curve is determined, the difference between each NDVI phenological curve and the average phenological curve can be calculated by Equation (4).
where n represents the number of points between the corresponding phenological phases of the reference curve and the target curve, and T represents the temporal resolution of the phenological curve.

Waveform Adjustment
Before extracting the differences between each NDVI phenological curve and the average phenological curve, it is necessary to preprocess the phenological curve of winter wheat, including waveform adjustment and polynomial fitting. Due to factors, such as temperature, latitude, and longitude, and the NDVI value of winter wheat, may be different in different regions even for the same phenological phase. For example, in some regions, when winter wheat enters the HD, the NDVI reaches 0.8, while in other regions, the NDVI only reaches 0.7. Therefore, it was necessary to adjust the NDVI time-series curve to reduce the influence of different amplitudes on the two phenological curves. Details of the waveform adjustment method are given in Appendix A.

Calculation of the Difference between NDVI Phenological Curves and the Average Phenological Curve
According to the time alignment relationship by the TWDTW algorithm in Section 3.1.2, when the phenological phase of the reference curve is known, the phenological phase of the target curve can be obtained. Figure 6 shows an example of how the TWDTW algorithm was used to determine the phenological phase. Two NDVI phenological curves were used, namely, a reference curve and a target curve, where the key phenological phases (GUD, HD, and MD) of the reference curve (termed GUDR, HDR, and MDR, respectively) are known. Three cases are analyzed here. In the first case, one point on the reference curve matches one point on the target curve; HDR and the location of HD on the target curve (HDT) belong to this situation; that is, there is only one line segment connecting HDR and HDT. This shows that, from the perspective of waveform similarity, HDT has similar characteristics to HDR, so HDT can be regarded as the HD of the target curve. In the second case, multiple points on the reference curve match one point on the target curve; GUDR and the location of GUD on the target curve (GUDT) belong to this situation. GUDR and the three points on the left match the GUDT point. From the perspective of alignment, although GUDT is aligned with four points, GUDR only matches GUDT, so GUDT can be regarded as the GUD of the target curve. In the third case, one point on the reference matches multiple points on the target curve; this is the case for MDR and the location of MD on the target curve (MDT). MDR is aligned with three points on the target curve. In this case, the average point is chosen as the MD of the target curve. After the phenological phase of the target curve is determined, the difference between each NDVI phenological curve and the average phenological curve can be calculated by Equation (4).
where n represents the number of points between the corresponding phenological phases of the reference curve and the target curve, and T represents the temporal resolution of the phenological curve.

Determination of GUD, HD, and MD
When the differences between each NDVI phenological curve and the average phenological curve (∆GUD, ∆HD, and ∆MD) were calculated, combined with the known average phenological phases. The phenological phases of the target curve were obtained according to Equations (5)- (7).
Using the same method, the GUD, HD, and MD were obtained for each winter-wheat NDVI phenological curve.

Accuracy Assessment
The coefficient of determination (R 2 ), root mean square error (RMSE), and bias were used to assess the accuracy of the values of GUD, HD, and MD. These metrics are defined as follows: where n represents the number of samples, y i represents the ground-observed phenological phase, y i represents the satellite-derived phenological phase, and y represents the average phenological phase.

Winter Wheat Distribution and Pure Winter-Wheat Pixels
The TWDTW method was used to extract winter wheat distribution in the Hebei area, and the result is shown in Figure 7. It can be seen from the figure that winter wheat in Hebei is mainly distributed in the central and southern plains of Hebei Province, and a small amount of cultivation is in the northeastern region of Hebei Province. Using 121 winter wheat and 107 non-winter wheat verification samples to verify the accuracy of the extraction results. The overall accuracy of the classification is 94.74%, and the kappa coefficient is 0.90. According to the winter-wheat area obtained by the TWDTW classification method, the area of winter wheat in Hebei Province in 2018 was found to be 2.5 million hectares, which is close to the value of 2.3 million hectares published in the Yearbook of Hebei Province 2018. Additionally, a total of 338 pure winter-wheat pixels were extracted (Figure 7). These pure winter-wheat pixels were used to generate the average phenological curve and analyze the spatial patterns of the phenological phases.

Spatial Patterns of GUD, HD, and MD
Based on the distribution of winter wheat, the TWDTW algorithm was used to extract three key phenological phases (GUD, HD, and MD) of winter wheat during the growth period in Hebei Province in 2018. The results for the GUD, HD, and MD are shown in Figure 8a-c, respectively, and the corresponding frequency distributions are shown in Figure 8d-

Spatial Patterns of GUD, HD, and MD
Based on the distribution of winter wheat, the TWDTW algorithm was used to extract three key phenological phases (GUD, HD, and MD) of winter wheat during the growth period in Hebei Province in 2018. The results for the GUD, HD, and MD are shown in Figure 8a, b, c, respectively, and the corresponding frequency distributions are shown in Figure  8d, e, and f, respectively. The GUD, HD, and MD all differ significantly between the south and north of Hebei Province.

Validation of the Extracted GUD, HD, and MD
As shown in Figure 9, the results of the phenological phase extraction were verified using the GUD, HD, and MD observed at 15 phenological monitoring stations. From Figure 9a, it can be seen that the phenological phase obtained by the TWDTW method is significantly correlated with the ground-observed phenological phase, with an R 2 greater than 0.59. The GUD has the highest R 2 of 0.71, while the MD has the lowest R 2 of 0.59. For the RMSE and bias, the HD has the smallest values of 5.72 days and 3.5 days, respectively, while the GUD has the largest values of 9.76 days and 8.4 days, respectively. As can be seen in Figure 9, the TWDTW method overestimates the GUD, HD, and MD. Accordingly, using the TWDTW method to extract the phenological phase has high accuracy and performs best for the HD.

Validation of the Extracted GUD, HD, and MD
As shown in Figure 9, the results of the phenological phase extraction were verified using the GUD, HD, and MD observed at 15 phenological monitoring stations. From Figure 9a, it can be seen that the phenological phase obtained by the TWDTW method is significantly correlated with the ground-observed phenological phase, with an R 2 greater than 0.59. The GUD has the highest R 2 of 0.71, while the MD has the lowest R 2 of 0.59. For the RMSE and bias, the HD has the smallest values of 5.72 days and 3.5 days, respectively, while the GUD has the largest values of 9.76 days and 8.4 days, respectively. As can be seen in Figure 9, the TWDTW method overestimates the GUD, HD, and MD. Accordingly, using the TWDTW method to extract the phenological phase has high accuracy and performs best for the HD.

Relationship between Winter-Wheat Phenology Pattern and Geographical Location
To investigate the distribution pattern of winter-wheat phenological phases, we analyzed the changes in the distribution of the main phenological phases of winter-wheat in Hebei Province with longitude and latitude. The results are shown in Figure 10. From this figure, it can be seen that the GUD, HD, and MD have a high correlation with latitude, with R 2 values of 0.658, 0.588, and 0.649, respectively. By establishing a linear regression equation, it was found that for each 1° increase in latitude, the GUD, HD, and MD are delayed by 4.84, 5.79, and 6.61 days, respectively (see Figure 10a, b, and c, respectively). The R 2 values of the GUD, HD, MD, with longitude are 0.258, 0.411, and 0.095, respectively (Figure 10d, e, and f, respectively). The HD has the strongest correlation with longitude. According to the regression equation, for each 1° increase in longitude, the HD is delayed by 5.2 days. The correlation between MD and longitude is the weakest, with an R 2 of only 0.095, indicating that MD does not change with longitude.

Relationship between Winter-Wheat Phenology Pattern and Geographical Location
To investigate the distribution pattern of winter-wheat phenological phases, we analyzed the changes in the distribution of the main phenological phases of winter-wheat in Hebei Province with longitude and latitude. The results are shown in Figure 10. From this figure, it can be seen that the GUD, HD, and MD have a high correlation with latitude, with R 2 values of 0.658, 0.588, and 0.649, respectively. By establishing a linear regression equation, it was found that for each 1 • increase in latitude, the GUD, HD, and MD are delayed by 4.84, 5.79, and 6.61 days, respectively (see Figure 10a-c, respectively). The R 2 values of the GUD, HD, MD, with longitude are 0.258, 0.411, and 0.095, respectively (Figure 10d-f, respectively). The HD has the strongest correlation with longitude. According to the regression equation, for each 1 • increase in longitude, the HD is delayed by 5.2 days. The correlation between MD and longitude is the weakest, with an R 2 of only 0.095, indicating that MD does not change with longitude.
To further analyze the influence of latitude and longitude on the phenological phases, the length of time between the GUD and HD (GUD-HD), between the HD and MD (HD-MD), and between the GUD and MD (GUD-MD) were extracted, and the relationships between these lengths and latitude and longitude were analyzed. The results are shown in Figure 11. As shown in the figure, the R 2 values between GUD-HD, HD-MD, and GUD-MD, and latitude are 0.071, 0.026, and 0.090, respectively (Figure 11a-c, respectively), which shows that the length of time between phenological phases varies little with latitude. Meanwhile, the R 2 values between GUD-HD, HD-MD, and GUD-MD, and longitude are also very small, namely, 0.213, 0.206, and 0.012, respectively (Figure 11d-f, respectively). This shows that the length of time between phenological phases varies little with longitude. The results show that the influence of longitude is greater than that of latitude on GUD-HD and HD-MD. In summary, changes in latitude and longitude impact the dates of the phenological phases, which is mainly due to the difference in temperature. As the latitude and longitude increase, the phenological phases are delayed, with latitude playing a larger role than longitude. However, latitude and longitude have little effect on the lengths of the phenological phases; across the study area, the lengths of time between GUD and HD, HD and MD, and GUD and MD are stable at 46, 41, and 87 days, respectively. longitude. The results show that the influence of longitude is greater than that of latitude on GUD-HD and HD-MD. In summary, changes in latitude and longitude impact the dates of the phenological phases, which is mainly due to the difference in temperature. As the latitude and longitude increase, the phenological phases are delayed, with latitude playing a larger role than longitude. However, latitude and longitude have little effect on the lengths of the phenological phases; across the study area, the lengths of time between GUD and HD, HD and MD, and GUD and MD are stable at 46, 41, and 87 days, respectively.   longitude. The results show that the influence of longitude is greater than that of latitude on GUD-HD and HD-MD. In summary, changes in latitude and longitude impact the dates of the phenological phases, which is mainly due to the difference in temperature. As the latitude and longitude increase, the phenological phases are delayed, with latitude playing a larger role than longitude. However, latitude and longitude have little effect on the lengths of the phenological phases; across the study area, the lengths of time between GUD and HD, HD and MD, and GUD and MD are stable at 46, 41, and 87 days, respectively.

Relationship between Winter-Wheat Phenology Pattern and Accumulated Temperature
In Hebei Province, winter wheat is generally sown in October, the GUD typically occurs at the end of February or early March of the following spring, the HD typically occurs at the end of March or the beginning of April, and the MD typically occurs at the end of May or early June. Therefore, the accumulated temperature was calculated from 01 October 2017 onwards: From 01 October 2017 to 28 February 2018 for the GUD; from 01 October 2017 to 30 April 2018 for the HD; and from 01 October 2017 to 31 May 2018 for the MD. The accumulated temperature between the GUD and HD is the accumulated temperature of the HD minus that of the GUD-that is, the accumulated temperature between 28 February and 30 April 2018. The accumulated temperature between the HD and MD was calculated between 30 April and 31 May 2018. The accumulated temperature between the GUD and MD was calculated between 28 February and 31 May 2018. By analyzing the relationships between the accumulated temperature and phenological phases, it was found that the correlations between the accumulated temperature and phenological phases are relatively high; the R 2 values between the GUD, HD, and MD, and the accumulated temperature in the corresponding period were 0.559, 0.506, and 0.595, respectively ( Figure 12a-c, respectively). The results showed that, with the decrease of accumulated temperature, the phenological phase is delayed. From Figure 12d-f, it can be seen that the correlation coefficients between the lengths of the phenological phases and the corresponding accumulated temperature are very low, with R 2 values of only 0.061, 0.033, and 0.072 for GUD, HD, and MD, respectively. perature of the HD minus that of the GUD-that is, the accumulated temperature between 28 February and 30 April 2018. The accumulated temperature between the HD and MD was calculated between 30 April and 31 May 2018. The accumulated temperature between the GUD and MD was calculated between 28 February and 31 May 2018. By analyzing the relationships between the accumulated temperature and phenological phases, it was found that the correlations between the accumulated temperature and phenological phases are relatively high; the R 2 values between the GUD, HD, and MD, and the accumulated temperature in the corresponding period were 0.559, 0.506, and 0.595, respectively (Figure 12a, b, and c, respectively). The results showed that, with the decrease of accumulated temperature, the phenological phase is delayed. From Figure 12d-f, it can be seen that the correlation coefficients between the lengths of the phenological phases and the corresponding accumulated temperature are very low, with R 2 values of only 0.061, 0.033, and 0.072 for GUD, HD, and MD, respectively.

Impact of the Calculation Method of the Average Phenological Curve on the Verification Results
Because the methods proposed in this paper involve using the average phenological curve and the target curve to calculate the difference between each NDVI phenological curve and the average phenological curve and extract the phenological phases, the average phenological curve will directly affect the accuracy of the results. Therefore, to analyze the influence of the calculation method for the average phenological curve on the accuracy of the phenological phase extraction, we calculated the average phenological curve using different randomly selected fractions of the 338 pure winter-wheat samples for sample

Impact of the Calculation Method of the Average Phenological Curve on the Verification Results
Because the methods proposed in this paper involve using the average phenological curve and the target curve to calculate the difference between each NDVI phenological curve and the average phenological curve and extract the phenological phases, the average phenological curve will directly affect the accuracy of the results. Therefore, to analyze the influence of the calculation method for the average phenological curve on the accuracy of the phenological phase extraction, we calculated the average phenological curve using different randomly selected fractions of the 338 pure winter-wheat samples for sample sizes between 5 and 100% with a 5% increment (i.e., 5%, 10%, 15%, . . . , 100%; that is, a total of 20 sample sizes) and used these average phenological curves to extract the GUD, HD, and MD using the TWDTW algorithm. The results are shown in Figure 13. As can be seen in the figure, changing the sample size had little effect on the R 2 , RMSE, and bias; for the GUD, HD, and MD, the R 2 was stable at about 0.7, 0.68, and 0.6, respectively; the RMSE remained stable at around 9.7, 5.7, and 7.0 days, respectively; and the bias was stable at around 8.4, 3.4, and 4.6 days, respectively.
Remote Sens. 2021, 13, x FOR PEER REVIEW 15 of 18 sizes between 5 and 100% with a 5% increment (i.e., 5%, 10%, 15%,…, 100%; that is, a total of 20 sample sizes) and used these average phenological curves to extract the GUD, HD, and MD using the TWDTW algorithm. The results are shown in Figure 13. As can be seen in the figure, changing the sample size had little effect on the R 2 , RMSE, and bias; for the GUD, HD, and MD, the R 2 was stable at about 0.7, 0.68, and 0.6, respectively; the RMSE remained stable at around 9.7, 5.7, and 7.0 days, respectively; and the bias was stable at around 8.4, 3.4, and 4.6 days, respectively. The above analysis suggests that using different numbers of winter-wheat samples to generate the average phenological curve has little effect on the results.

Comparison with Other Methods
To further verify the performance and advantages of using the TWDTW algorithm to extract winter-wheat phenology, several other methods were also used to extract the GUD, HD, and MD of winter wheat. These methods were the RT10 and RT20 thresholding Figure 13. The values of (a) the coefficient of determination (R 2 ), (b) root mean square error (RMSE), and (c) bias obtained for average phenological curves produced using different sample sizes.
The above analysis suggests that using different numbers of winter-wheat samples to generate the average phenological curve has little effect on the results.

Comparison with Other Methods
To further verify the performance and advantages of using the TWDTW algorithm to extract winter-wheat phenology, several other methods were also used to extract the GUD, HD, and MD of winter wheat. These methods were the RT10 and RT20 thresholding methods and a method involving the maximum of the curvature change rate of the fitted logistic curve (CCR max ) [13]; RT10 and RT20 are typical dynamic thresholding methods and RT10 and RT20 is a change detection method. The method of extracting HD uses the maximum method [35]. The method of extracting MD uses the RT10 method [24].
We compared the results of the phenological extraction using the TWDTW algorithm with those using the RT10, RT20, and CCR max methods, as illustrated in Figure  14. Compared with the other methods, in terms of R 2 , the TWDTW algorithm does not offer a notable advantage (Figure 14a). For the GUD, the TWDTW algorithm is slightly better than the RT20 and CCR MAX methods, and slightly inferior to the RT10 method. For the HD, the performance of the TWDTW algorithm is similar to that of the CCR MAX method. For the MD, the TWDTW algorithm performs better than the RT10 method. However, from the perspective of the RMSE and bias (Figure 14b,c, respectively), for the GUD, the performance of the TWDTW algorithm is similar to those of the RT10 and CCR MAX methods and significantly better than that of the RT20 method. That is, for the HD and MD, the TWDTW algorithm is significantly better than the CCR MAX and RT10 methods, respectively. Figure 13. The values of (a) the coefficient of determination (R 2 ), (b) root mean square error (RMSE), and (c) bias obtained for average phenological curves produced using different sample sizes.

Comparison with Other Methods
To further verify the performance and advantages of using the TWDTW algorithm to extract winter-wheat phenology, several other methods were also used to extract the GUD, HD, and MD of winter wheat. These methods were the RT10 and RT20 thresholding methods and a method involving the maximum of the curvature change rate of the fitted logistic curve (CCRmax) [13]; RT10 and RT20 are typical dynamic thresholding methods and RT10 and RT20 is a change detection method. The method of extracting HD uses the maximum method [35]. The method of extracting MD uses the RT10 method [24].
We compared the results of the phenological extraction using the TWDTW algorithm with those using the RT10, RT20, and CCRmax methods, as illustrated in Figure 14. Compared with the other methods, in terms of R 2 , the TWDTW algorithm does not offer a notable advantage (Figure14a). For the GUD, the TWDTW algorithm is slightly better than the RT20 and CCRMAX methods, and slightly inferior to the RT10 method. For the HD, the performance of the TWDTW algorithm is similar to that of the CCRMAX method. For the MD, the TWDTW algorithm performs better than the RT10 method. However, from the perspective of the RMSE and bias (Figure14b and c, respectively), for the GUD, the performance of the TWDTW algorithm is similar to those of the RT10 and CCRMAX methods and significantly better than that of the RT20 method. That is, for the HD and MD, the TWDTW algorithm is significantly better than the CCRMAX and RT10 methods, respectively.

Conclusions
In this paper, we proposed a novel phenological phase determination method based on the TWDTW algorithm using MODIS NDVI time-series data. The advantages of the TWDTW algorithm are as follows: (1) Based on extracted pure winter-wheat pixels, an average phenological curve is generated, which can effectively reduce the influence of mixed pixels; (2) the algorithm completes the time alignment of the entire growth period and can extract GUD, HD, and MD simultaneously. The phenological phase extraction using the TWDTW algorithm was found to have high accuracy. In terms of verification accuracy, the RMSEs of the GUD, HD, and MD were 9.76, 5.72, and 6.98 days, respectively. Additionally, the proposed TWDTW method was also demonstrated to perform better than several other methods. Furthermore, it was shown that, in Hebei Province, the GUD, HD, and MD are mainly affected by latitude and accumulated temperature; as the latitude increases from south to north, the GUD, HD, and MD are delayed. For each 1 • rise in latitude, the GUD, HD, and MD are delayed by 4.84, 5.79, and 6.61 days, respectively. The higher the accumulated temperature, the earlier the phenological phases occurred. However, the latitude and accumulated temperature had little effect on the length of the phenological phases. Throughout Hebei Province, the lengths of time between GUD and HD, HD and MD, and GUD and MD were stable at 46, 41, and 87 days, respectively. The S B_NEW (t)= ⎩ ⎨ ⎧ S B1 (t)* B max -B min1 +A max -B max * B max -B min1 32 ≤ t ≤ T A S B2 (t)* A max -A min2 B max -B min2 +A max -B max * A max -A min2 B max -B min2 T A < t ≤ 186 (13) where Bmax represents the maximum of SB(t), and Bmin1 and Bmin2 represent the minimum of SB1(t) and SB2(t), respectively. Amax represents the maximum value of SA(t), and Amin1 and Amin2 represent the minimum value of SA1(t)and SA2(t), respectively. Figure A1. A comparison of the NDVI time-series curve and the average phenological curve before (a) and after (b) waveform adjustment. Figure A1. A comparison of the NDVI time-series curve and the average phenological curve before (a) and after (b) waveform adjustment.