Mapping Winter Wheat in North China Using Sentinel 2A / B Data: A Method Based on Phenology-Time Weighted Dynamic Time Warping

: Accurate mapping of winter wheat over a large area is of great signiﬁcance for guiding policy formulation related to food security, farmland management, and the international food trade. Due to the complex phenological features of winter wheat, the cloud contamination in time-series imagery, and the inﬂuence of the soil / snow background on vegetation indices, there remains no e ﬀ ective method for mapping winter wheat at a medium spatial resolution (10–30 m). In this study, we proposed a novel method called phenology-time weighted dynamic time warping (PT-DTW) for identifying winter wheat based on Sentinel 2A / B time-series data. The main advantages of PT-DTW include (1) the use of phenological features in two periods, i.e., the greenness increase before winter and greenness decrease after heading, which are common to all winter wheat and are distinct from the features of other land cover types, and (2) the use of the normalized di ﬀ erential phenology index (NDPI) instead of traditional vegetation indices to provide more robust vegetation information and to suppress the adverse impacts of soil and snow cover, especially during the before-winter growth period. The proposed PT-DTW method was employed for winter wheat mapping based on Sentinel 2A / B data on the Huang-Huai Plain, China. Validation with visually interpreted samples showed that the produced winter wheat map achieved an overall classiﬁcation accuracy of 89.98% and a kappa coe ﬃ cient of 0.7978, outperforming previous winter wheat classiﬁcation methods. Moreover, the planting area derived from PT-DTW agreed well with census data at the municipal level, with a coe ﬃ cient of determination of 0.8638, indicating that the winter wheat map produced at 20 m resolution was reliable overall. Therefore, the PT-DTW method is recommended for winter wheat mapping over large areas.


Introduction
Wheat is a grain crop that is cultivated worldwide and provides nearly 20% of all calories consumed due to its strong adaptability to various temperature and water conditions [1]. As the country with the greatest wheat consumption worldwide, China maintains a vast wheat cultivation area to ensure its food security [2]. The winter wheat cultivation region in North China is the largest wheat-producing area in China and is experiencing considerable effects from climate change, wheat-producing area in China and is experiencing considerable effects from climate change, alterations in food preferences, and the international food trade [3,4]. Accurate and timely mapping of winter wheat in the area is indispensable for formulating policies with respect to food security, agriculture structure adjustment, and the international food trade. Remote sensing techniques therefore play a unique role in monitoring the planting areas of winter wheat for intensive observation in both spatial and temporal respects.
Winter wheat has unique phenological features, including several important phenology stages during its life cycle: sowing, seedling, tillering, overwintering, greening-up, jointing, heading, grouting, maturity, and harvesting. Such a growth process can be observed by the curve of the vegetation index (VI) time series (called the VI curve below) derived from remotely sensed data ( Figure 1a). Generally, the VI curve of winter wheat exhibits a unique pattern with two peaks, at the tillering and heading stages. Most remote sensing-based mapping methods use this phenological feature [5][6][7][8][9][10][11][12] and have achieved satisfactory results in small areas. However, the unique phenological features of winter wheat can be interfered and become inconspicuous or obscured in larger areas with various agroclimatic conditions. For example, the VI curves display large shifts and distortions across large areas (Figure 1b as an example). In addition, the VI peak at the tillering stage before winter is usually not obvious due to the varying growth status after tillering and the effects of the soil background and snow cover. Accordingly, accurate winter wheat mapping over a large area is still a major challenge due to the high intraclass variability of winter wheat [13][14][15][16][17].  [18] curves and the corresponding phenological stages of winter wheat. (b) NDPI curves of winter wheat from different locations in North China. NDPI curves in (a,b) are derived from sentinel-2 MSI data. Curve D is from a location of a rotation of winter wheat and summer corn.
Generally, remote sensing-based methods for winter wheat mapping can be grouped into three types: classification-based methods, index-based methods, and curve similarity-based methods. The first type of method is the classification of one or two satellite images acquired at particular phenology stages, such as at two peak occurrence periods [5][6][7][8][9][10][11]. Unfortunately, the opportunity to acquire high-quality images covering a large area within a specific period is usually limited by long satellite revisit cycles and cloud contamination. Moreover, the collection of a large number of training data covering a large area not only is time consuming, but also has a high field cost. Accordingly, winter wheat mapping over large areas using classification-based methods is not a reasonable solution.
In the second type of method, the VI curve is employed to design specific indices for winter wheat based on its unique phenological features. For example, Pan et al. [8] developed a crop proportion phenology index (CPPI) that used enhanced vegetation index (EVI) values related to the seedling, tillering, heading, and harvesting stages. The CPPI was demonstrated to have a good ability to estimate the winter wheat fraction in a 250 m MODIS pixel through experiments in two small agricultural regions in Beijing and Jiangsu provinces of China. Tao et al. [19] used mathematical transformation to generate the peak before winter feature (PBWF) index associated with the EVI peak amplitude before winter to map winter wheat on the North China Plain, although PBWF neglects the phenological feature variability over large areas. Qiu et al. [14] also proposed two winter wheat indices representing ascending and descending trends before and after the heading date, while ignoring the VI increase before winter. These winter wheat indices typically require prior knowledge of important phenology dates for each pixel due to the large phenological shifts over large areas (see Figure 1b). To address this issue, based on phenological observation data from agrometeorological stations, Qiu et al. [14] built a regression model for the heading date (also the early growing period length) with latitude and altitude, with which the important phenology dates in a specific location can be estimated. Nevertheless, the accuracy of the regression model remains inadequate because the spatial variation in phenology dates is also affected by factors other than latitude and altitude. Meanwhile, such prior knowledge collection or estimation is usually difficult without intensive field surveys or auxiliary data support, which hinders the wide application of this kind of method over large areas [13,15,17].
In the third group of methods, winter wheat is identified by a VI curve similarity measurement without requiring any prior knowledge of important phenology stages. Sun et al. [9] used the Euclidean distance to measure the similarity between any EVI2 curve and a reference EVI2 curve of winter wheat, although the Euclidean distance has limited tolerance for the shifts and distortion of winter wheat curves. Zhang et al. [10] mapped winter wheat in Luoyang city, China with a Kullback-Leibler divergence (KLD)-based method, which was proven effective for land cover classification [20], whereas the KLD classifier lacks verification of winter wheat mapping over larger areas.
Recently, the dynamic time warping (DTW) algorithm, which was originally proposed for speech recognition [21], has also been used in land cover classification [22][23][24][25] due to its ability to tolerate both the shifts and distortions of VI curves. The original DTW identifies the best alignment between two VI curves based on the Euclidean distance without any constraints, which often leads to unreasonable alignment. To address this issue, warping time constraints and a distortion penalty are often imposed to restrict the warping range. For example, Petitjean et al. [26] introduced a constant elapsed time constraint into DTW. Maus et al. [27] proposed time-weighted DTW (TWDTW) with a penalty strategy. The TWDTW method has been successfully utilized for multicrop mapping [28]. Due to the sufficient flexibility of DTW and the considerable efforts devoted to improving DTW performance, the DTW algorithm has been successfully employed for land cover and crop type mapping over large areas [27][28][29][30][31]. Unfortunately, DTW has been applied to winter wheat mapping in only a very limited number of studies [28], possibly due to the more complex characteristics of winter wheat VI curves. The application of DTW to winter wheat mapping over large areas faces the following challenges: (1) how to utilize the special phenological features of winter wheat and (2) how to improve the quality of VI data, especially in the early development stages during which VIs are easily affected by soil background and snow cover. To address these two challenges, a new curve matching method was developed in this study, termed phenology-time weighted dynamic time warping (PT-DTW), which specifically considers the phenological features of winter wheat. In addition, the normalized differential phenology index (NDPI), also called NDVI+ [18,32], was employed instead of a traditional vegetation index (e.g., NDVI or EVI) to mitigate the effects of the snow and soil background. With the proposed method, winter wheat on the Huang-Huai Plain in North China was mapped at 20 m resolution based on Sentinel 2A/B data, considering the high resolution of Sentinel 2A/B data in both the spatial and temporal dimensions.

Study Area and Datasets
The Huang-Huai Plain, located at 29.08 • N-42.67 • N latitude and 105.18 • E-122.70 • E longitude (Figure 2), is the largest wheat production zone in China [33] and accounts for 44% of the wheat planting area and 60% of the total wheat production in the country [34]. In this region, there is commonly a rotation between winter wheat and corn/paddy rice, which means that winter wheat is sown after corn/paddy rice harvesting and that corn/paddy rice is planted after the harvesting of winter wheat [2]. Winter wheat is sown in September or October and harvested in May or June of the following year [5]. As the climatic conditions vary widely over the study area, shown by the 2015 annual average air temperature (http://www.resdc.cn/data.aspx?DATAID=228) in Figure 2, there is large intraclass variability in the VI curves of winter wheat in different locations ( Figure 1b). In general, the seedling date advances and the heading date is delayed from south to north, leading to a large distortion in the VI curves of winter wheat across the study area. In addition, from north to south, the before-winter peak becomes less obvious due to the higher temperatures in winter in the southern region. Sentinel 2 is a European wide-swath, high-resolution, multispectral imaging mission. Two identical satellites (Sentinel 2A/B) provide remotely sensed data with a short revisit cycle (2-3 days at mid-latitudes) and high spatial resolution (10, 20, and 60 m for different spectral bands). A total of 184 tiles of Sentinel 2A/B images cover the entire study area, and all images with less than 80% cloud coverage from September 2017 to June 2018 were downloaded (https://earthexplorer.usgs.gov/). All of the images were used to reconstruct high-quality NDPI curves throughout the growth period of winter wheat (see Section 2.2.1). Sentinel 2 is a European wide-swath, high-resolution, multispectral imaging mission. Two identical satellites (Sentinel 2A/B) provide remotely sensed data with a short revisit cycle (2-3 days at Remote Sens. 2020, 12, 1274 5 of 22 mid-latitudes) and high spatial resolution (10, 20, and 60 m for different spectral bands). A total of 184 tiles of Sentinel 2A/B images cover the entire study area, and all images with less than 80% cloud coverage from September 2017 to June 2018 were downloaded (https://earthexplorer.usgs.gov/). All of the images were used to reconstruct high-quality NDPI curves throughout the growth period of winter wheat (see Section 2.2.1).
The ground truth data were interpreted through high-resolution images acquired at 495 sites in five subregions for winter wheat mapping and accuracy assessment ( Figure 3). The layout of these sites was designated along north-south and east-west transects to ensure coverage of multiple agroclimatic zones. The high-resolution images were acquired by an unmanned aerial vehicle system (DJI FC6310) in the growing season of 2017-2018 or Google Earth imagery captured in the winter of 2017. At each sampling site, the corresponding winter wheat patches were visually interpreted using the high-resolution imagery. Then, the interpreted maps were projected to the corresponding 20 m Sentinel pixels ( Figure 3). In total, 62,239 ground-truth pixels with 20 m spatial resolution were collected for accuracy evaluation including 27,669 winter wheat and 34,570 other pixels (including forest, barren land, artificial surface, vegetable, paddy rice, water, and grass/shrub land). For the determination of method parameters 80 samples from pure pixels were selected (Section 2.2.3.1), and all the others were used for accuracy assessment. In addition, agricultural census data from 61 municipalities of Beijing, Henan, Anhui, Shandong, Shanxi, Shaanxi, and Anhui in 2017 were acquired from the EPS China Data Platform (http://www.epschinadata.com/). The census data were used to evaluate the overall reliability of the identification of planting areas of winter wheat at the municipal level.

Methods
Based on Sentinel 2A/B NDPI time series, we proposed a PT-DTW method for winter wheat mapping in the study area. The workflow consists of the following steps: (a) reconstructing high-quality NDPI time series; (b) optimizing the parameters of the PT-DTW method using simulated samples; and (c) identifying winter wheat and performing accuracy assessment ( Figure 4).

Methods
Based on Sentinel 2A/B NDPI time series, we proposed a PT-DTW method for winter wheat mapping in the study area. The workflow consists of the following steps: (a) reconstructing high-quality NDPI time series; (b) optimizing the parameters of the PT-DTW method using simulated samples; and (c) identifying winter wheat and performing accuracy assessment ( Figure 4).

Reconstruction of High-quality NDPI Curves
NDPI curves were derived from the multispectral bands of Sentinel 2A/B MSI data in this step. Firstly, the Sentinel 2 data were atmospherically corrected by Sen2Cor (http://step.esa.int/main/third-party-plugins-2/sen2cor/). Secondly, the images of the red (band 4), near-infrared (band 8), and shortwave infrared (band11) bands were all resampled to 20 m spatial resolution, the channel information was listed in Table 1. The NDPI was thus calculated as follows: where α was set to 0.74, its most effective value for suppressing the variability of the soil and snow backgrounds [18,32]. Thirdly, the NDPI curve was composited with a seven-day maximum composing criterion. We used all of the images with less than 80% cloud contamination, but there was still a lack of observations on some dates and in some areas. To address this issue, the missing observations were linearly interpolated from the valid values on the two nearest dates before and after the missing observations. Finally, Savitzky-Golay filtering for vegetation index reconstruction [35] was performed to reduce the atmospheric effects further and to generate the seven-day NDPI curves at 20 m spatial resolution.

Reconstruction of High-quality NDPI Curves
NDPI curves were derived from the multispectral bands of Sentinel 2A/B MSI data in this step. Firstly, the Sentinel 2 data were atmospherically corrected by Sen2Cor (http://step.esa.int/main/ third-party-plugins-2/sen2cor/). Secondly, the images of the red (band 4), near-infrared (band 8), and shortwave infrared (band11) bands were all resampled to 20 m spatial resolution, the channel information was listed in Table 1. The NDPI was thus calculated as follows: where α was set to 0.74, its most effective value for suppressing the variability of the soil and snow backgrounds [18,32]. Thirdly, the NDPI curve was composited with a seven-day maximum composing criterion. We used all of the images with less than 80% cloud contamination, but there was still a lack of observations on some dates and in some areas. To address this issue, the missing observations were linearly interpolated from the valid values on the two nearest dates before and after the missing observations. Finally, Savitzky-Golay filtering for vegetation index reconstruction [35] was performed to reduce the atmospheric effects further and to generate the seven-day NDPI curves at 20 m spatial resolution.

Development of the PT-DTW Method
As mentioned before, the phenological features of winter wheat varied widely over large areas ( Figure 1b). More specifically, winter wheat is sown later and harvested earlier and thus experiences a shorter growing period from north to south because of the warmer climatic conditions in the south. There is additional and obvious intraclass variation induced by the different growing statuses in winter. In colder areas, winter wheat stops growing and goes dormant in winter and becomes green in early spring. This pattern creates a U-shaped valley in the NDPI curve. Whereas in warmer areas, winter wheat grows continuously through winter and early spring, displaying an insignificant valley in the NDPI curve. These complex features of the NDPI curves may explain why DTW-based methods are not widely used for winter wheat mapping. Fortunately, the winter wheat across the study area also shares common NDPI characteristics in two periods (Figure 5a), i.e., increase before winter and decrease after heading (see the areas with the light yellow background in Figure 5a). These features are not only common features of all winter wheat, but are also distinct from the features of other land cover types. Accordingly, it is still possible to apply the DTW-based method in winter wheat mapping over larger areas if these common phenological features can be enhanced. Therefore, the proposed PT-DTW imposes greater weights on these two important periods when calculating the DTW distance ( Figure 5b). The details of the PT-DTW are as follows. As mentioned before, the phenological features of winter wheat varied widely over large areas ( Figure 1b). More specifically, winter wheat is sown later and harvested earlier and thus experiences a shorter growing period from north to south because of the warmer climatic conditions in the south. There is additional and obvious intraclass variation induced by the different growing statuses in winter. In colder areas, winter wheat stops growing and goes dormant in winter and becomes green in early spring. This pattern creates a U-shaped valley in the NDPI curve. Whereas in warmer areas, winter wheat grows continuously through winter and early spring, displaying an insignificant valley in the NDPI curve. These complex features of the NDPI curves may explain why DTW-based methods are not widely used for winter wheat mapping. Fortunately, the winter wheat across the study area also shares common NDPI characteristics in two periods (Figure 5a), i.e., increase before winter and decrease after heading (see the areas with the light yellow background in Figure 5a). These features are not only common features of all winter wheat, but are also distinct from the features of other land cover types. Accordingly, it is still possible to apply the DTW-based method in winter wheat mapping over larger areas if these common phenological features can be enhanced. Therefore, the proposed PT-DTW imposes greater weights on these two important periods when calculating the DTW distance ( Figure 5b). The details of the PT-DTW are as follows.

Classical DTW Algorithm
The classical DTW algorithm aims to search for optimal alignment of the target and reference curves, which correspond to the warping with the shortest distance between the two curves. Suppose that there are two curves,

Classical DTW Algorithm
The classical DTW algorithm aims to search for optimal alignment of the target and reference curves, which correspond to the warping with the shortest distance between the two curves. Suppose that there are two curves, U = (u 1 , u 2 , u 3 , . . . , um) and R = (r 1 , r 2 , r 3 , . . . , r n ), denoting the target and reference curves (Figure 6a), respectively. An m×n matrix dm × n ( Figure 6a) is constructed to denote the distances between the pairs of elements in the two curves: where each element in the matrix dm × n is defined as d i,j = ui − r j . Then, an accumulated distance matrix D is computed by a recursive sum of the minimal distances, such that Remote Sens. 2020, 12, 1274 8 of 22 with the following initial conditions:

Logistic Time Penalty
Considering that an excessively large distortion is unreasonable for NDPI curve matching, logistic TWDTW [27] was utilized to avoid the excessively large distortions in classical DTW. A strategy was used that assigns a penalty weight to each element in the distance matrix: where Mi,j is the penalty weight following a logistic function with midpoint β and steepness α : We set α = 0.1 and β = 100 according to the recommendation of Maus et al. [27]. With this penalty on the distance, the unreasonably large distortion of classic DTW can be avoided.

PT-DTW
PT-DTW is specifically designed for winter wheat identification by increasing the importance of the two feature periods and decreasing the impacts of other periods. Specifically, the distances associated with the two feature periods (VI increase before winter and VI decrease after heading) ( Figure 5b) were given greater weights, and the other distances were given smaller weights. The phenological weighted function was designed as shown in Equation (9) where N1 is the total number of alignment associated with the two feature periods of the reference curve; N2 is the total number of alignment associated with the nonfeature periods; ω is the parameter controlling the weight allocation, which was optimized in the simulation experiment (Section 2.2.3); j is the time-phase index in the reference NDPI curve; and Phases 8, 16, 33, and 41 correspond to the Then, the warping path P = (P 1 , P 2 , . . . , P L ) in matrix d is determined through the reverse algorithm ( Figure 6a). P preserves the indices of the matching elements between the two curves, also called alignment of two curves ( Figure 6b). P starts with P 1 = (1, 1) and ends with P L = (m, n). The searching algorithm of P is as follows: Finally, the shortest warping path P is acquired, and the corresponding DTW distance between U and R is: where L is the length of P. Equation (3) guarantees the monotonicity and continuity condition [29].
Equations (4) and (6) guarantee that the shortest warping path is restricted to starting with 'beginning to beginning' and terminating with 'ending to ending'.

Logistic Time Penalty
Considering that an excessively large distortion is unreasonable for NDPI curve matching, logistic TWDTW [27] was utilized to avoid the excessively large distortions in classical DTW. A strategy was used that assigns a penalty weight to each element in the distance matrix: Remote Sens. 2020, 12, 1274 where M i,j is the penalty weight following a logistic function with midpoint β and steepness α: We set α = 0.1 and β = 100 according to the recommendation of Maus et al. [27]. With this penalty on the distance, the unreasonably large distortion of classic DTW can be avoided.

PT-DTW
PT-DTW is specifically designed for winter wheat identification by increasing the importance of the two feature periods and decreasing the impacts of other periods. Specifically, the distances associated with the two feature periods (VI increase before winter and VI decrease after heading) (Figure 5b) were given greater weights, and the other distances were given smaller weights. The phenological weighted function was designed as shown in Equation (9): 16] or j ∈ [33,41] (1 − ω)/N2, j ∈ other bands (9) where N 1 is the total number of alignment associated with the two feature periods of the reference curve; N 2 is the total number of alignment associated with the nonfeature periods; ω is the parameter controlling the weight allocation, which was optimized in the simulation experiment (Section 2. A smaller PT-DTW distance indicates a higher similarity between the reference and target curve, and winter wheat pixels can be distinguished from others if the distance is less than a certain threshold (T; Equation (11)): The threshold was also optimized in the simulation experiment (Section 2.2.3).

Determination of PT-DTW Parameters
Several parameters need to be determined for the use of PT-DTW classification, including the NDPI reference curve of winter wheat, phenological weight (ω), and threshold (T). Typically, parameter optimization relies on large training samples. To avoid a large sample collection effort, we designed a simulation experiment to expand the small sample set to a large simulated dataset by using a linear spectral mixing (LSM) model [36,37]. As shown in Figure 4b, 80 manually interpreted samples of different land cover types were firstly collected, one of which was determined to be the reference NDPI curve; then, a large simulated sample set was generated by using a linear mixture model; finally, the parameters ω and T were optimized such that they yielded the highest classification accuracy for the simulated dataset. The details are as follows.

Collection of Samples and Determination of the NDPI Reference Curve for Winter Wheat
Based on prior knowledge of the study area, eight land cover types were considered: winter wheat, paddy rice, vegetable, water, artificial surface, forest, grass/shrub land, and barren land. Ten NDPI curves for each land cover type were selected across the study area with reference to high-resolution Google Earth imagery. In total, 80 NDPI curves were selected as initial samples Regarding the NDPI reference curve, an ideal reference curve should lie in the center of all winter wheat curves. Therefore, the NDPI reference curve for winter wheat was determined as the one with the minimum average Euclidean distance from the nine other winter wheat curves (Figure 7). model; finally, the parameters and T were optimized such that they yielded the highest classification accuracy for the simulated dataset. The details are as follows.

Collection of Samples and Determination of the NDPI Reference Curve for Winter Wheat
Based on prior knowledge of the study area, eight land cover types were considered: winter wheat, paddy rice, vegetable, water, artificial surface, forest, grass/shrub land, and barren land. Ten NDPI curves for each land cover type were selected across the study area with reference to high-resolution Google Earth imagery. In total, 80 NDPI curves were selected as initial samples Regarding the NDPI reference curve, an ideal reference curve should lie in the center of all winter wheat curves. Therefore, the NDPI reference curve for winter wheat was determined as the one with the minimum average Euclidean distance from the nine other winter wheat curves ( Figure  7).

Sample Set Expansion by Using an LSM Model
Considering that the initial sample set, which included only 80 NDPI curves (Figure 8a), was too small to optimize the PT-DTW parameters effectively, the sample set was expanded to a large dataset by using a linear spectral mixing (LSM) model ( Figure 8). Mixed NDPI curves were simulated by mixing the NDPI curves of winter wheat endmembers with the NDPI curves of other land cover endmembers. Here, a three-endmember mixture model was considered because it covered most mixing scenarios.
NDPI mix = f wheat NDPI wheat + f nonwheat1 NDPI nonwheat1 + f nonwheat2 NDPI nonwheat2 subject to f wheat + f nonwheat1 + f nonwheat2 = 1 (12) where f wheat is the fraction of wheat endmember, and f nonwheat1 and f nonwheat2 are the fractions of the two nonwheat endmembers. f wheat , f nonwheat1 , and f nonwheat2 were generated from the probability density functions shown in Figure 8b. Such distributions guarantee a uniform distribution (0-1) of the wheat fraction and sum-to-one constraint for the three endmembers. Details about the generation process are introduced in the supplementary material. In total, 20,000 mixed NDPI curves were simulated, among which 10,000 were positive samples (f wheat > 50%) and the other 10,000 were negative samples (f wheat < 50%) as shown in Figure 8c.

Parameter Optimization
We varied ω from 0.1 to 1.0 in increments of 0.1 and varied T from the minimum value to the maximum value in 300 steps. Under each ω, the minimum and maximum T are the median PT-DTW distances of the positive and negative samples, respectively. Then, we classified the simulated dataset by using PT-DTW and calculated the overall accuracy with each combination of ω and T (Figure 9). With this optimization procedure, the optimal values of ω = 1 and T = 0.0765, corresponding to the highest overall accuracy (OA; 92.85%), were acquired. As shown in Figure 9b, a large class separability between the positive and negative samples was acquired when ω was equal to 1, indicating that it is a good proxy for distinguishing winter wheat from other land cover types. nonwheat endmembers. fwheat, fnonwheat1, and fnonwheat2 were generated from the probability density functions shown in Figure 8b. Such distributions guarantee a uniform distribution (0-1) of the wheat fraction and sum-to-one constraint for the three endmembers. Details about the generation process are introduced in the supplementary material. In total, 20,000 mixed NDPI curves were simulated, among which 10,000 were positive samples (fwheat > 50%) and the other 10,000 were negative samples (fwheat < 50%) as shown in Figure 8c.

Parameter Optimization
We varied from 0.1 to 1.0 in increments of 0.1 and varied T from the minimum value to the maximum value in 300 steps. Under each , the minimum and maximum T are the median PT-DTW distances of the positive and negative samples, respectively. Then, we classified the simulated dataset by using PT-DTW and calculated the overall accuracy with each combination of and T (Figure 9). With this optimization procedure, the optimal values of = 1 and T = 0.0765, corresponding to the highest overall accuracy (OA; 92.85%), were acquired. As shown in Figure 9b, a large class separability between the positive and negative samples was acquired when was equal to 1, indicating that it is a good proxy for distinguishing winter wheat from other land cover types.

Winter Wheat Mapping and Accuracy Assessment
With the parameters determined above, the PT-DTW distance to the NDPI reference curve of winter wheat was calculated. Finally, the winter wheat pixels in the study area were extracted with the optimal threshold (T = 0.0765). The TWDTW method [27], CBAH method [14], and KLD method [10] were used to map the winter wheat for comparison. The thresholds in the TWDTW, CBAH, and

Winter Wheat Mapping and Accuracy Assessment
With the parameters determined above, the PT-DTW distance to the NDPI reference curve of winter wheat was calculated. Finally, the winter wheat pixels in the study area were extracted with the optimal threshold (T = 0.0765). The TWDTW method [27], CBAH method [14], and KLD method [10] were used to map the winter wheat for comparison. The thresholds in the TWDTW, CBAH, and KLD methods were also optimized by performing similar simulation experiments. The overall accuracies and kappa coefficients of these four methods were calculated based on the visually interpreted samples. In addition, the winter wheat areas extracted using the four methods were compared with agricultural census data at the municipal level.

Sentinel 2-Derived Winter Wheat Map for 2017-2018 by PT-DTW Method
Based on Sentinel 2A/B data, a winter wheat map at 20 m resolution was produced using the proposed PT-DTW method. As shown in Figure 10, winter wheat was mainly distributed in Henan and Jiangsu Provinces, Northern Anhui Province, Southern Hebei Province, and Western Shandong Province, which is generally consistent with the results of previous studies [14,18]. The zoomed-in images in Figure 10 show rich spatial details of the winter wheat cropland, and the cultivated parcels are easily recognizable. Considering the fragmented crop cultivation in North China, a 20-m-resolution winter wheat map can provide valuable information for cropland management and policy making.

Quantitative Evaluation
The classification accuracies of the PT-DTW, TWDTW, CBAH, and KLD methods were evaluated with 62,239 ground-truth pixels, including 27,669 winter wheat and 34,570 other pixels. Table 2 shows the confusion matrixes and classification accuracy indices of these four methods. The proposed PT-DTW method provides the highest accuracy. The CBAH method also achieves satisfactory accuracy, whereas the TWDTW and KLD methods yield poor classification accuracy.

Quantitative Evaluation
The classification accuracies of the PT-DTW, TWDTW, CBAH, and KLD methods were evaluated with 62,239 ground-truth pixels, including 27,669 winter wheat and 34,570 other pixels. Table 2 shows the confusion matrixes and classification accuracy indices of these four methods. The proposed PT-DTW method provides the highest accuracy. The CBAH method also achieves satisfactory accuracy, whereas the TWDTW and KLD methods yield poor classification accuracy. The winter wheat pixels derived from Sentinel 2A/B data were aggregated to the municipal level for comparison with the census data, i.e., the municipal-level seeded area of winter wheat in 2017. The cultivated areas of winter wheat estimated using the PT-DTW method matched best with the agricultural census data at the municipal level with a coefficient of determination of 0.8638 and a mean absolute error of 858.8412 km 2 , although slight overestimation can be found in the municipalities with larger areas of winter wheat (Figure 11a). In contrast, the TWDTW method shows significant overestimation in the municipalities with small areas of winter wheat, and the KLD method shows the largest deviation from the census data. There is also obvious overestimation of the areas estimated by CBAH, especially in the municipalities with small areas of winter wheat, despite good correlations with the census data in other municipalities.

Comparison of Spatial Details
Typical mapping examples of four locations were compared to show the differences among the PT-DTW, CBAH, KLD, and TWDTW methods ( Figure 12). In general, the PT-DTW-derived maps show good consistency with the distribution of winter wheat, whereas the other methods produced more classification errors to different degrees. As shown in Figure 12I, II, the KLD-derived winter wheat map is weakly consistent with the spatial pattern of winter wheat patches and has large omission and commission errors. Two shortcomings resulted in the poor performance of the KLD method in winter wheat mapping. Firstly, without reasonable alignment between the target and reference NDPI curves, the KL divergence could not suppress the winter wheat variability across large areas. Secondly, the KLD is a general measure of curve similarity and does not include any consideration of the specific phenology features of winter wheat. The CBAH method is much better for winter wheat mapping. The unique phenological characteristics of winter wheat (i.e., the VI ascending and descending trends before and after heading) are utilized in CBAH classification. However, CBAH overlooks the greenness increase before winter and misclassifies some spring-sown and summer-harvested crops as winter wheat (Figure 12 III). Moreover, latitude and altitude are not sufficient to capture all of the variations in sowing and heading date, which are also probably influenced by the local climate and field management. As an example, in Figure 12 IV, the predicted heading date (DOY 2018113) is earlier than the true heading date (DOY 2018136), leading to some omission error in winter wheat mapping. The TWDTW method, without the consideration of specific phenology periods of winter wheat, yielded poor classification accuracy. TWDTW misclassified forests in mountainous areas (Figure 12 II) and spring-sown crops (Figure 12 III) as winter wheat and produced salt and pepper noise in the patch of winter wheat cropland (Figure 12 I). Compared to the published methods, PT-DTW produced the best winter wheat map due to its effective alignment and its consideration of the specific phenology periods in curve matching.

Comparison of Spatial Details
Typical mapping examples of four locations were compared to show the differences among the PT-DTW, CBAH, KLD, and TWDTW methods ( Figure 12). In general, the PT-DTW-derived maps show good consistency with the distribution of winter wheat, whereas the other methods produced more classification errors to different degrees. As shown in Figure 12I, II, the KLD-derived winter

Superiority of NDPI
The NDPI is a newly developed VI that has the advantages of strengthening vegetation information and suppressing soil color variation and snow cover compared to traditional Vis [18,32]. Early growth before winter is the important phenology feature used in PT-DTW to identify winter wheat. However, as winter wheat grows weakly with a low leaf area during that period, traditional VIs may be influenced considerably by the snow and soil background. Therefore, the NDPI was used instead of the NDVI and EVI to capture the weak greenness signal of winter wheat more effectively. To verify whether the NDPI outperforms other VIs in capturing the increase in greenness before winter, the NDPI curves were compared with the NDVI, EVI, and EVI2 curves based on the visually interpreted winter wheat samples. As shown in Figure 13a, the average NDPI

Superiority of NDPI
The NDPI is a newly developed VI that has the advantages of strengthening vegetation information and suppressing soil color variation and snow cover compared to traditional Vis [18,32]. Early growth before winter is the important phenology feature used in PT-DTW to identify winter wheat. However, as winter wheat grows weakly with a low leaf area during that period, traditional VIs may be influenced considerably by the snow and soil background. Therefore, the NDPI was used instead of the NDVI and EVI to capture the weak greenness signal of winter wheat more effectively. To verify whether the NDPI outperforms other VIs in capturing the increase in greenness before winter, the NDPI curves were compared with the NDVI, EVI, and EVI2 curves based on the visually interpreted winter wheat samples. As shown in Figure 13a, the average NDPI curve of winter wheat shows a more obvious increase before winter than the curve of the other VIs. To express the magnitude of the increase quantitatively, a before-winter growth index (BWGI) was calculated for the different VIs: (13) where VI max is the maximum VI of winter wheat during the winter growing season and VI min is the minimum VI before the winter peak. Specifically, we identified the maximum VI between 01 November 2017 and 31 December 2017 as VI max and the minimum VI between 01 October 2017 and 31 November 2017 as VI min . As shown in Figure 13b, the BWGI of NDPI is significantly larger than those of NDVI, EVI, and EVI2 (p < 0.01 in a t-test), indicating that NDPI captured the winter wheat growth before winter more effectively than NDVI, EVI, and EVI2.
Remote Sens. 2020, 03, x FOR PEER REVIEW 18 of 25 curve of winter wheat shows a more obvious increase before winter than the curve of the other VIs.
To express the magnitude of the increase quantitatively, a before-winter growth index (BWGI) was calculated for the different VIs: max min BWGI VI VI = − (13) where max VI is the maximum VI of winter wheat during the winter growing season and min VI is the minimum VI before the winter peak. Specifically, we identified the maximum VI between 01 November 2017 and 31 December 2017 as max VI and the minimum VI between 01 October 2017 and 31 November 2017 as min VI . As shown in Figure 13b, the BWGI of NDPI is significantly larger than those of NDVI, EVI, and EVI2 (p < 0.01 in a t-test), indicating that NDPI captured the winter wheat growth before winter more effectively than NDVI, EVI, and EVI2.

Benefits of the PT-DTW Method
To illustrate the benefits of PT-DTW quantitatively, we further compared the two-class separability, as well as the interclass and intraclass variability, of the PT-DTW, TWDTW, and traditional Euclidean distances of the simulated and visually interpreted samples. All three distances measure the similarity between two curves based on Euclidean distances, while TWDTW and PT-DTW improve the Euclidean distance by including specific curve alignment and phenological feature weighting considerations. To verify the effectiveness of the two abovementioned enhancements, the two-class separability is defined as follows [38,39]: where w μ and o μ are the mean values for the winter wheat and other samples, respectively; wheat and other land cover types, denoting the intraclass variability. As shown in Figure 14, the PT-DTW distance produces the most distinctive histograms between winter wheat and other land cover types among the three distances, followed by the TWDTW and Euclidean distances. Table 3 provides more detailed information about the two-class separability, interclass variability, and intraclass variability. The TWDTW and PT-DTW distances show much smaller intraclass variability

Benefits of the PT-DTW Method
To illustrate the benefits of PT-DTW quantitatively, we further compared the two-class separability, as well as the interclass and intraclass variability, of the PT-DTW, TWDTW, and traditional Euclidean distances of the simulated and visually interpreted samples. All three distances measure the similarity between two curves based on Euclidean distances, while TWDTW and PT-DTW improve the Euclidean distance by including specific curve alignment and phenological feature weighting considerations. To verify the effectiveness of the two abovementioned enhancements, the two-class separability is defined as follows [38,39]: where µ w and µ o are the mean values for the winter wheat and other samples, respectively; µ w − µ o denotes the interclass variability; and σ w and σ o are the standard deviations of winter wheat and other land cover types, denoting the intraclass variability. As shown in Figure 14, the PT-DTW distance produces the most distinctive histograms between winter wheat and other land cover types among the three distances, followed by the TWDTW and Euclidean distances. Table 3 provides more detailed information about the two-class separability, interclass variability, and intraclass variability.
The TWDTW and PT-DTW distances show much smaller intraclass variability of winter wheat than the Euclidean distance (Table 3), indicating that DTW-based alignment can effectively tolerate the shifts and distortions of the NDPI curves of winter wheat. Despite its effectiveness in suppressing the intraclass variability of winter wheat, the classification accuracy of TWDTW is only slightly higher than that of the Euclidean distance, because the interclass variability is also reduced. With the optimization of phenology weighting, PT-DTW greatly increased the interclass variability while simultaneously maintaining low intraclass variability ( Table 3), suggesting that weighting on certain phenology periods successfully enhanced the distinctive features of winter wheat. In summary, the above results confirmed the effectiveness of TWDTW alignment and phenology weighting in extracting winter wheat with PT-DTW.
Remote Sens. 2020, 03, x FOR PEER REVIEW 19 of 25 of winter wheat than the Euclidean distance ( Table 3), indicating that DTW-based alignment can effectively tolerate the shifts and distortions of the NDPI curves of winter wheat. Despite its effectiveness in suppressing the intraclass variability of winter wheat, the classification accuracy of TWDTW is only slightly higher than that of the Euclidean distance, because the interclass variability is also reduced. With the optimization of phenology weighting, PT-DTW greatly increased the interclass variability while simultaneously maintaining low intraclass variability ( Table 3), suggesting that weighting on certain phenology periods successfully enhanced the distinctive features of winter wheat. In summary, the above results confirmed the effectiveness of TWDTW alignment and phenology weighting in extracting winter wheat with PT-DTW.

Transferability of PT-DTW Method
Winter wheat phenology is dominated by climatic factors and field management. Thus, the phenology shifts appear not only over large areas, but also across different years. To validate whether the PT-DTW method can be applied to other years, one tile of the Sentinel 2 (50SLH) NDPI time series from 2018-2019 was used to derive the winter wheat map through PT-DTW with the same parameters employed in the 2017-2018 mapping. The winter wheat cultivation in this area located in Hebei province, experienced some changes because of the new crop rotation and fallow policy implemented since 2016 [40,41]. As shown in Figure 15, the transitions between the winter wheat and fallow land were successfully detected by the PT-DTW method without retraining the PT-DTW parameters. The land cover types of 3786 cloud-free pixels were interpreted visually through images acquired in 2017/12/24 and 2018/12/09, considering that winter wheat shows

Transferability of PT-DTW Method
Winter wheat phenology is dominated by climatic factors and field management. Thus, the phenology shifts appear not only over large areas, but also across different years. To validate whether the PT-DTW method can be applied to other years, one tile of the Sentinel 2 (50SLH) NDPI time series from 2018-2019 was used to derive the winter wheat map through PT-DTW with the same parameters employed in the 2017-2018 mapping. The winter wheat cultivation in this area located in Hebei province, experienced some changes because of the new crop rotation and fallow policy implemented since 2016 [40,41]. As shown in Figure 15, the transitions between the winter wheat and fallow land were successfully detected by the PT-DTW method without retraining the PT-DTW parameters. The land cover types of 3786 cloud-free pixels were interpreted visually through images acquired in 2017/12/24 and 2018/12/09, considering that winter wheat shows greenness within cropland in winter. The change detection accuracy was then evaluated, as shown in Table 4. An OA of 93.66% in change detection indicates the great potential of PT-DTW for rapid winter wheat mapping in multiple years. Although some of the published methods, such as CBAH, can achieve comparable accuracy with the support of auxiliary data, as mentioned in Sections 3.2 and 3.3, the PT-DTW method is proven to be more flexible and practical, given the lower demand for auxiliary data and great temporal transferability.
Remote Sens. 2020, 03, x FOR PEER REVIEW 20 of 25 greenness within cropland in winter. The change detection accuracy was then evaluated, as shown in Table 4. An OA of 93.66% in change detection indicates the great potential of PT-DTW for rapid winter wheat mapping in multiple years. Although some of the published methods, such as CBAH, can achieve comparable accuracy with the support of auxiliary data, as mentioned in Sections 3.2 and 3.3, the PT-DTW method is proven to be more flexible and practical, given the lower demand for auxiliary data and great temporal transferability. In addition to the temporal transferability of winter wheat mapping, the proposed PT-DTW has the potential to be modified for mapping other crops besides winter wheat. Although PT-DTW cannot be directly used to map other crops because the parameters were specifically determined based on the phenological features of winter wheat, it would not be difficult to set up a new group of parameters to classify other crops using a strategy similar to that employed in this study. If key phenological stages could be determined based on prior knowledge, the parameters could also be trained through a small number of field samples using the method described in Section 2.2.3. After the parameters were redetermined, the retrained PT-DTW could be used to classify other crops.  In addition to the temporal transferability of winter wheat mapping, the proposed PT-DTW has the potential to be modified for mapping other crops besides winter wheat. Although PT-DTW cannot be directly used to map other crops because the parameters were specifically determined based on the phenological features of winter wheat, it would not be difficult to set up a new group of parameters to classify other crops using a strategy similar to that employed in this study. If key phenological stages could be determined based on prior knowledge, the parameters could also be trained through a small number of field samples using the method described in Section 2.2.3. After the parameters were redetermined, the retrained PT-DTW could be used to classify other crops.

Influence of Cloud Contamination
Cloud contamination is a considerable obstacle for time-series matching. Thus, a simulation experiment was conducted to explore the sensitivity of the PT-DTW method to cloud contamination.
First, all the visually interpreted samples with little cloud contamination were selected as the benchmark data for analysis. Second, cloudy NDPI curves were simulated by randomly replacing the clear observations with cloudy ones (NDPI was set to 0). The different cloud percentages and different fractions of cloud observations in key phenological stages in the total cloud observations were considered (details in the Supplementary Material). Finally, the simulated cloudy NDPI time-series were used to classify winter wheat, by PT-DTW with the same time series reconstruction procedure (Section 2.2.1) and parameters. As shown in Figure 16a, the OA decreases marginally when the cloud percentage is less than 50%, indicating that the Savizky-Golay filter performed well in reconstructing the NDPI time series. Figure 16b shows that the uneven cloud distributions in key and nonkey phenological stages also marginally affect the classification accuracy when the cloud percentage is less than 50%, which may be because the overall NDPI curve shape did not change considerably regardless of how the cloud observations were distributed in key and nonkey phenological stages. The above results indicate that the PT-DTW method is suitable for winter wheat mapping when the percentage of clear Sentinel 2 observations is greater than 50% throughout the growing period of winter wheat.

Limitations of the PT-DTW Method
Several limitations remain in the proposed PT-DTW method. Firstly, the Euclidean distance is calculated to measure the difference between the target and reference curves, leading to a relatively high sensitivity to the magnitude variability of NDPI curves. Thus, winter wheat pixels with poor growth status are likely to be misclassified due to their lower NDPI values. Induction of the distances into non-Euclidean space is a possible solution; however, this approach would increase the mathematical complexity of the algorithms. The data quality may also affect the mapping accuracy. As mentioned in the previous section, PT-DTW performs poorly when the cloud percentage of Sentinel 2 data is greater than 50%, which could occur in Southern China. The fusion of cloud-free SAR data could provide VI time-series data with better quality and help improve the performance of wheat mapping [42].

Conclusions
In this study, a PT-DTW method was proposed for winter wheat mapping over a large area based on the NDPI curves derived from Sentinel 2A/B data. With specific improvements over the

Limitations of the PT-DTW Method
Several limitations remain in the proposed PT-DTW method. Firstly, the Euclidean distance is calculated to measure the difference between the target and reference curves, leading to a relatively high sensitivity to the magnitude variability of NDPI curves. Thus, winter wheat pixels with poor growth status are likely to be misclassified due to their lower NDPI values. Induction of the distances into non-Euclidean space is a possible solution; however, this approach would increase the mathematical complexity of the algorithms. The data quality may also affect the mapping accuracy. As mentioned in the previous section, PT-DTW performs poorly when the cloud percentage of Sentinel 2 data is greater than 50%, which could occur in Southern China. The fusion of cloud-free SAR data could provide VI time-series data with better quality and help improve the performance of wheat mapping [42].

Conclusions
In this study, a PT-DTW method was proposed for winter wheat mapping over a large area based on the NDPI curves derived from Sentinel 2A/B data. With specific improvements over the classical DTW and TWDTW, PT-DTW can tolerate the variability of VI curves of winter wheat over large areas and works well without considerable sampling effort. Comparison experiments showed that the PT-DTW method outperforms previous methods, achieves reasonable mapping accuracy, and is transferable to winter wheat mapping in different years. Therefore, the PT-DTW method is recommended for winter wheat mapping over large areas and in multiple years, requiring only a small number of samples for parameter training.
With the proposed PT-DTW, the winter wheat of North China in 2017-2018 was mapped at 20 m resolution, which is the highest resolution in a winter wheat map of this area to our knowledge. The map produced at this spatial resolution can provide rich spatial details for the area, which has fragmented crop cultivation.