Comparison of Remote Sensing Time-Series Smoothing Methods for Grassland Spring Phenology Extraction on the Qinghai–Tibetan Plateau

: Accurate evaluation of start of season (SOS) changes is essential to assess the ecosystem’s response to climate change. Smoothing method is an understudied factor that can lead to great uncertainties in SOS extraction, and the applicable situation for di ﬀ erent smoothing methods and the impact of smoothing parameters on SOS extraction accuracy are of critical importance to be clariﬁed. In this paper, we use MOD13Q1 normalized di ﬀ erence vegetation index (NDVI) data and SOS observations from eight agrometeorological stations on the Qinghai–Tibetan Plateau (QTP) during 2001–2011 to compare the SOS extraction accuracies of six popular smoothing methods (Changing Weight (CW), Savitzky-Golay (SG), Asymmetric Gaussian (AG), Double-logistic (DL), Whittaker Smoother (WS) and Harmonic Analysis of NDVI Time-Series (HANTS)) for two types of di ﬀ erent SOS extraction methods (dynamic threshold (DT) with 9 di ﬀ erent thresholds and double logistic (Zhang)). Furthermore, a parameter sensitivity analysis for each smoothing method is performed to quantify the impacts of smoothing parameters on SOS extraction. Finally, the suggested smoothing methods and reference ranges for the parameters of di ﬀ erent smoothing methods were given for grassland phenology extraction on the QTP. The main conclusions are as follows: (1) the smoothing methods and SOS extraction methods jointly determine the SOS extraction accuracy, and a bad denoising performance of smoothing method does not necessarily lead to a low SOS extraction accuracy; (2) the default parameters for most smoothing methods can result in acceptable SOS extraction accuracies, but for some smoothing methods (e.g., WS) a parameter optimization is necessary, and the optimal parameters of the smoothing method can increase the R 2 and reduce the RMSE of SOS extraction by up to 25% and 331%; (3) The main inﬂuencing factor of the SOS extraction using the DT method is the stability of the minimum value in the NDVI curve, and for the Zhang method the curve shape before the peak of the NDVI curve impacts the most; (4) HANTS is the most stable method no matter with (ﬁtness = 35.05) or without parameter optimization (ﬁtness = 33.52), which is recommended for QTP grassland SOS extraction. The ﬁndings of this study imply that remote sensing-based vegetation phenology extraction can be highly uncertain, and a careful selection and parameterization of the time-series smoothing method should be taken to achieve an accurate result.


Introduction
As a result of anthropogenic activities, global warming has brought widespread impacts on the terrestrial ecosystem [1,2]. Vegetation phenology is a sensitive indicator of climate change, and it has been widely reported that global warming has altered the vegetation phenology in the past few decades [3][4][5]. Changes in vegetation phenology can bring changes to the interaction between the biosphere and the atmosphere, and finally lead to changes in carbon balance [6], water balance [7] and even vegetation feedback mechanisms towards climate change [8]. Therefore, an accurate assessment of vegetation phenology is of crucial importance for understanding the impact of global climate change on terrestrial ecosystem cycle of carbon and for making effective adaptive management decisions [9].
There are two main data sources of vegetation phenology observations, one is from ground and near-ground measurements, and the other one is generated from remote sensing data [10][11][12]. Traditional ground and near-ground measurements have the advantage of high accuracy, but they lack the spatial continuity due to the limited number of observation stations [5]. Remote sensing, on the contrary, can provide simultaneous and spatial-continuous observations at large scales, and thus can present more insights at a spatial scale [13]. In recent years, satellite-based remote sensing technology has been widely used for vegetation phenology extraction and monitoring [9,[14][15][16][17][18][19][20]. However, the extraction of vegetation phenology using remote sensing technology can have serious uncertainties due to different data sources, time-series smoothing methods and phenology extraction methods. These uncertainties may even lead to different conclusions regarding the same question. For instance, using GIMMS NDVI data, Yu et al. [21] found an advanced followed by a delayed trend for start of season (SOS) on the Qinghai-Tibetan Plateau (QTP). However, Zhang et al. [22] reported that the SOS of alpine vegetation on the QTP is continuously advanced from 1982 to 2011. Only a few studies have been conducted to quantify the uncertainties regarding remote sensing phenology extraction, and these studies mainly focus on the raw remote sensing data and phenology extraction methods [23][24][25][26][27][28]. However, we have not yet acquired an explicit understanding of the impact of remote sensing time-series smoothing methods as well as the parameters of the smoothing methods on SOS extraction.
The current remote sensing time-series smoothing methods have large differences in the model structure [29], which may lead to great differences among smoothed curves and, furthermore, the extracted vegetation phenology. At present, the research regarding smoothing methods mostly focuses on the quality of curve reconstruction rather than the extraction of vegetation phenology [30,31].
Although several studies compared the phenology extraction results based on different smoothing methods [32,33], there is still a lack in quantitative evaluation of the interaction of different smoothing methods and different phenology extraction methods for grassland phenology extraction. Another problem regarding the time-series smoothing for phenology extraction is the parameter setting of the smoothing methods. Most smoothing methods require setting smoothing parameters manually [34][35][36], although the default parameter values were suggested for some of the smoothing methods when they were proposed (e.g., the default values of parameters m and d for Savitzky-Golay filter were suggested to be 4 and 6, respectively [35]), the optimal parameter values can vary across different study area and vegetation types due to different vegetation growth trajectories, and improper smoothing parameter values will lead to great uncertainties for the smoothing results [37][38][39] and inaccurate phenology extraction results. The vegetation phenology on the QTP has been widely studied due to its nature of being vulnerable and sensitive to climate change [21,24,40]. However, to our knowledge, there is still no comprehensive research of smoothing parameter comparison for grassland phenology extraction on the QTP; therefore, it is difficult for the users to parameterize the smoothing methods without an expert knowledge of this region.
In this study, we focus on quantifying the uncertainties arise from remote sensing time-series smoothing methods as well as the parameters of the smoothing methods for SOS extraction. To solve this problem, we used MODIS normalized difference vegetation index (NDVI) data and ground phenology measurement records from eight agrometeorological stations to compare six smoothing Remote Sens. 2020, 12, 3383 3 of 26 methods for grassland SOS extraction. The Qinghai-Tibetan Plateau, one of the most sensitive areas to global climate change, was chosen as the study area to conduct the study. This study aimed to examine three questions: (1) how do smoothing methods and their parameters impact on phenology extraction? (2) What are the main factors that bring uncertainties to the phenology extraction using remote sensing data? (3) What are the applicable conditions for different smoothing methods?

Study Area
The Qinghai-Tibetan Plateau stands in the center of Eurasia, spanning 31 longitude degrees (73 • 18 52"E-104 • 46 59"E) and 13 latitude degrees (26 • 00 12"N-39 • 46 50"N), with a total area of almost 2.5 × 10 6 km 2 [41]. The QTP is known as the roof of the world and the third pole (the highest plateau in the world), with an average altitude of 4000 m. As for the climate, the QTP is the coldest region at the same latitude due to its high altitude, the average surface temperature of the coldest month is between −10 to −15 • C. The average daily radiation on the QTP is 21 MJ m −2 day −1 , which is much higher than that of other areas at the same altitude [42]. The vegetation type on the QTP varies according to different temperature and water conditions, presenting a landscape of forests-grassland-desert from southeast to northwest ( Figure 1). The Alpine steppe has the widest distribution and the largest area among all vegetation types, covering 59.15% of the total area of the QTP [40]. The cold, dry and strong radiation conditions make the QTP one of the most sensitive areas to global climate change [43], and the QTP has experienced significant warming with the surface temperature increased by~1.8 • C during the past five decades [44], making it a hotspot for studying phenology change and its relationship with climate change [24,40,45].
Remote Sens. 2018, 5, x FOR PEER REVIEW  3 of 26 to global climate change, was chosen as the study area to conduct the study. This study aimed to examine three questions: (1) how do smoothing methods and their parameters impact on phenology extraction? (2) What are the main factors that bring uncertainties to the phenology extraction using remote sensing data? (3) What are the applicable conditions for different smoothing methods?

Study Area
The Qinghai-Tibetan Plateau stands in the center of Eurasia, spanning 31 longitude degrees (73°18′52″ E-104°46′59″ E) and 13 latitude degrees (26°00′12″ N-39°46′50″ N), with a total area of almost 2.5 × 10 6 km 2 [41]. The QTP is known as the roof of the world and the third pole (the highest plateau in the world), with an average altitude of 4000 m. As for the climate, the QTP is the coldest region at the same latitude due to its high altitude, the average surface temperature of the coldest month is between −10 to −15 °C. The average daily radiation on the QTP is 21 MJ m −2 day −1 , which is much higher than that of other areas at the same altitude [42]. The vegetation type on the QTP varies according to different temperature and water conditions, presenting a landscape of forestsgrassland-desert from southeast to northwest ( Figure 1). The Alpine steppe has the widest distribution and the largest area among all vegetation types, covering 59.15% of the total area of the QTP [40]. The cold, dry and strong radiation conditions make the QTP one of the most sensitive areas to global climate change [43], and the QTP has experienced significant warming with the surface temperature increased by ~1.8 °C during the past five decades [44], making it a hotspot for studying phenology change and its relationship with climate change [24,40,45].

Ground Observation Data
The grassland phenology records for the QTP were collected from the nationwide phenological observation network established by the China Meteorological Administration, and the dataset involves observations from eight agrometeorological stations. We excluded data that are not herbal, and the final dataset includes seven phenological metrics (Green-up date, Tillering data, Heading date, Florescence, Senescence) for nine herbaceous plants (Elymus nutans Griseb., Poa pratensis, Stipa krylovii Roshev, Leymus secalinus (Georgi) Tzvel., Agropyron cristatum, Festuca, Kobresia humilis (C. A. Mey ex Trauvt.) Sergievskaya, Astragalus adsurgen, Artemisia scoparia). We used the green-up date to represent the SOS date, and if a station recorded more than 1 herbal species then the averaged SOS was used. The spatial distribution of these eight stations is shown in Figure 1.

Remote Sensing Data and Processing
In this study, we used MODerate resolution Imaging Spectroradiometer (MODIS) Vegetation Indices product (MOD13Q1) version 6 during 2001-2011 for remote sensing-based SOS extraction, MOD13Q1 provides composited vegetation index (NDVI and enhanced vegetation index) time-series data at a 16-day interval and 250 m resolution, and we used the NDVI data for its wide-usage for SOS extraction [46][47][48] and for that most time-series smoothing methods are proposed based on NDVI data [35,36,49]. MODIS Land Cover Type product (MCD12Q1) version 6 was used for the extraction of grassland pixels, and the 500 m land cover data was resampled into 250 m resolution using nearest neighboring to consist with the NDVI data.
We used the grassland NDVI pixels within the 10 km × 10 km (43 pixel × 43 pixel at 250 m resolution) spatial range around each agrometeorological station to represent the grassland for the agrometeorological stations. Firstly, MOD13Q1 and MCD12Q1 data within the 10 km × 10 km window for each station were downloaded from the MODIS/VIIRS Global Subsets (https://modis.ornl.gov/cgibin/MODIS/global/subset.pl). Secondly, the MOD13Q1 NDVI time-series data within the 10 km × 10 km window for each station were smoothed, and the SOS for each pixel were extracted. Finally, the grassland pixels in the window were identified according to MCD12Q1 land cover data, and the averaged SOS of the grassland pixels in the window is used as the final SOS extracted using remote sensing data of the corresponding station.

The Time-Series Smoothing Methods
We chose six popular smoothing methods for the comparison based on previous studies [25,[30][31][32]50,51], including one empirical method (Changing-weight (CW)), four curve fitting methods (Savitzky-Golay (SG), Asymmetric Gaussian (AG) and Double-logistic (DL) and Whittaker Smoother (WS)) and one data transformation methods (Harmonic Analysis of NDVI Time-Series (HANTS)). Although the parameters of smoothing methods varied for different applications, since there is no similar work, we test the accuracy of default values (derived from the original or high citation papers) of each method. The details for each smoothing method and default values are shown in Table 1.

The Phenology Extraction Methods
We used two widely-adopted vegetation phenology extraction methods for grassland phenology extraction, one is the vegetation index ratio threshold method (also known as the dynamic threshold method), which has received a broad attention and application [23,55,56]; the other phenology extraction method is the double-logistic method (also known as the Zhang method), which has been reported to have ecologically meaning [15] and also has been applied in many phenology extraction related works [18,57,58].

Dynamic Threshold Method
The dynamic threshold method was originally proposed by White et al. [59]. This method normalizes the NDVI time-series data of a single pixel to 0-1, and uses the percentage of NDVI to represent the vegetation growth status of the pixel range.
where NDVI t is the NDVI value of t-th year and NDVI max is the maximum NDVI value of NDVI t , while the NDVI min is the minimum of the left half curve for SOS extraction. When NDVI ratio exceeds a certain threshold, the corresponding day of year (DOY) is determined as the SOS. Different thresholds represent different phenology phases. The thresholds chosen by previous studies are mostly 20% and 50% [21,26]. In this paper, we used thresholds ranging from 10% to 50% with a 5% increment, and there are 9 thresholds in total.

Double-Logistic Method
Double-logistic method is proposed by Zhang et al. [15], and hence the method is also called the Zhang method. In this method, the logistic curve is used to simulate the vegetation growth curve by piecewise fitting. After simulating the vegetation growth using the logistic model, the curvature change rate of the fitted logistic model is used to determine the vegetation phenology. The model's rate of change of curvature (RCC) is: where z = ea + bz, a and b are the fitting parameters, d is the initial background VI value, and c + d is the maximum VI value. SOS is the date in which the first local maximum of the curvature change rate of the model.

The Evaluation of the Phenology Extraction Accuracy
We evaluated the accuracy of the phenology extraction by comparing the satellite-derived SOS against the ground observed SOS. A simple linear regression was conducted where the ground observed SOS was used as the independent variable and the satellite-derived SOS was used as the predictor. The coefficient of determination (R 2 ) and the root mean square error (RMSE) of the linear regression were calculated to represent the trend similarity and the numerical closeness between the satellite-derived SOS and the ground observed SOS, respectively. Moreover, we defined a metric "fitness" to comprehensively measure the accuracy of the phenology extraction by combining R 2 and RMSE, which is calculated using the following equation: A lower fitness illustrates a smaller phenology extraction error, or a higher phenology extraction accuracy. Based on the temporal interval of MOD13Q1 (16 days), a fitness equal or less than 64 (corresponding to R 2 ≥ 0.25 [60] and RMSE ≤ 16 [21]) is defined as a good result.

Sensitivity Analysis of the Smoothing Parameters
Different parameter values directly affect the smoothing performance and the final SOS extraction accuracy. In this study, we used the grid search method to test the SOS extraction accuracy under different smoothing parameters. Each parameter set for different smoothing methods was tested at a stepwise according to the parameter ranges and steps in Table 1, and the SOS was extracted based on the corresponding parameter combination, and a fitness was used as the measure of the SOS extraction accuracy. The grid search was conducted for the 4 parameterized smoothing methods (CW, SG, WS and HANTS) combined with the two types of phenology extraction method (9 thresholds of DT and the Zhang method). For each phenology extraction method, the smoothing methods with the parameter sets that achieved the lowest fitness were chosen as the optimal smoothing methods, and the four optimal methods are abbreviated as O-CW, O-SG, O-WS, and O-HANTS.
To further quantify the contribution of each smoothing parameter to the SOS extraction accuracy, a standardized multi-linear regression was conducted by taking the value of different smoothing parameters as independent variables and the fitness as the dependent variable. The contribution of each parameter to the SOS extraction accuracy was characterized with the absolute value of the corresponding coefficient.

Statistical Significance Test of the Results
In order to investigate whether the value of SOS extracted by different smoothing methods for the same extraction method is significantly different, and whether the accuracy of the same smoothing method for different extraction methods is significantly different, we performed a two-way ANOVA analysis and least significant difference (LSD) post hoc multiple comparisons using the smoothing method and extraction method as the factors, and SOS as the dependent.

Performance of the Time-Series Smoothing Results
The main differences among smoothing methods are mainly reflected in the smoothing performance towards the low value, the ability of keeping the maximum value of the curve and the shape of the fitted curve. For default parameters of the smoothing methods (including AG and DL), the smoothing methods managed to eliminate most of the noise in the raw NDVI curves, but different smoothing methods resulted in smoothed curves with varied maximum and minimum values. In all the smoothing methods, SG always managed to maintain the annual maximum value, while the annual maximum value of WS is significantly smaller than that of the other methods ( Overall, the smoothing methods result in smoothed curved with similar shapes and trends under the condition with almost no noise. However, as the noise increases, the shapes of the curves smoothed by different methods starts to differ, and the curves smoothed by WS has the greatest difference compared with other methods, especially in the first half of the curve, with lower starting values and larger slopes before DOY 150. It can be easily observed that HANTS is sensitive to small fluctuations in the curve during the smoothing process, resulting in the presence of multiple peaks in the NDVI curve. Such peaks usually appear at the beginning of the NDVI curve, but the heights of these peaks are relatively low (Figure 3a,c). For the case with continuous noise, most methods show good resistance against the continuous noise, but NDVI curves smoothed by SG and CW show much lower value than the raw NDVI right after the appearance of the continuous noise ( Figure 3c). from the 9 DTs for further analysis. It can be found that although the tendency variation of the minimum values among the 4 smoothing methods (AG and DL excluded) (average SD = 0.18 for DT and 0.16 for Zhang method) are still much greater than that of the maximum values (average SD = 0.04 for DT and 0.03 for Zhang method) after adopting the optimal parameters, the tendency variation of the minimum values decreases substantially compared with that with the default parameters (decrease by 51.4% for DT and 56.76% for Zhang method), and there are much fewer points with extremely low values and fewer fluctuations of the minimum values in the smoothed curves ( Figures  2a and 4a,b). O-HANTS have the smoothest effects for low values in the raw NDVI curve, resulting in an annual minimum NDVI value (0.23) higher than that of other methods.  As for the optimal parameters, the data range of the parameters is similar for the DT and Zhang method (Tables 2 and 3). Since we focus on the optimal parameter of the smoothing methods that achieve the highest accuracy, we only use the DT with the lowest fitness for each smoothing method from the 9 DTs for further analysis. It can be found that although the tendency variation of the minimum values among the 4 smoothing methods (AG and DL excluded) (average SD = 0.18 for DT and 0.16 for Zhang method) are still much greater than that of the maximum values (average SD = 0.04 for DT and 0.03 for Zhang method) after adopting the optimal parameters, the tendency variation of the minimum values decreases substantially compared with that with the default parameters (decrease by 51.4% for DT and 56.76% for Zhang method), and there are much fewer points with extremely low values and fewer fluctuations of the minimum values in the smoothed curves (Figures 2a and 4a,b). O-HANTS have the smoothest effects for low values in the raw NDVI curve, resulting in an annual minimum NDVI value (0.23) higher than that of other methods.

Inter-Comparison of SOS Extraction Using Different Smoothing Methods
The default parameters of all the smoothing methods result in an overestimation of SOS regarding DT in station 1, 2 and 7, which is especially severe for station 1 and 2, with an overestimation of SOS of ~20 days (Figure 2b). However, for Zhang method, almost all the smoothing methods with the default parameter result in prediction of SOS close to the observation level except for WS, which shows a pervasive underestimation of SOS of 35.79 ± 16.51 for all the stations and years (Figure 2c). It is interesting to note that there are four points (marked with red circles in Figure 2b,c) with large underestimation (up to ~80 days) of SOS for both DT and Zhang SOS extraction methods, and most of the smoothing methods show this underestimation with varying degrees, in which AG and DL have the greatest underestimations. As for the SOS extraction accuracy, most of the smoothing methods achieve a RMSE less than 16 days for both DT and Zhang methods (Tables 2 and  3), which is acceptable according to the temporal interval of the NDVI data [21], but WS with default parameters results in a RMSE that is even larger than 2 times the temporal interval of the NDVI data ( Table 3). It is worth noting that even with the default parameters, all the smoothing methods have different performances for different SOS extraction methods, and for most of the smoothing methods, the fitness for Zhang is significantly lower than that for DT, with an average fitness 75.36 lower for Zhang compared with that for DT. For the two phenology extraction methods, the smoothing methods that lead to the highest (HANTS) and lowest (WS) accuracy are the same, in which WS has much greater RMSE for Zhang method, while the R 2 of HANTS SOS almost doubled for Zhang method compared with DT.   The optimal parameters improve the SOS extraction accuracies for all the smoothing methods compared with the default parameters, and the fitness decrease from 84. 38 (Tables 2  and 3), and the average fitness of the four optimal smoothing methods for DT and Zhang are 89.04 (decreased by 43.24% compared with the default parameters) and 41.68 (decreased by 77.49% compared with the default parameters), respectively. WS has the greatest improvement in SOS extraction accuracy for both DT and Zhang after adopting the optimal parameters, with a RMSE decrease and R 2 increase of 27.81 (3.82) and 24% (19%) for Zhang (DT), respectively. Moreover, the underestimation of SOS of WS for Zhang method has been significantly improved after using the optimal parameters (Figures 2c and 4b).
For the DT methods with nine thresholds for SOS extraction, O-CW, AG and DL result in similar fitness that are above 150, while the fitness for O-SG are generally lower than the 3 methods, but the values are still all over 100 ( Figure 5). On the contrary, O-WS and O-HANTS has much smaller fitness compared with the other 4 smoothing methods for the 9 DTs, with an average fitness value less than 90 for the 9 DTs (Table A1). The low fitness of O-HANTS and O-WS can be attributed to both higher R 2 and smaller RMSE compared with other smoothing methods (average R 2 and RMSE are 0.23 (0.23) and 13.75 (19.30) for O-WS (O-HANTS)). However, O-HANTS has higher R 2 compared with O-WS for the best threshold, while O-WS shows small RMSE (RMSE < 12) with a wide range of thresholds (15%-40%), leading to a close fitness for these thresholds. The best thresholds of DT for grassland SOS ranges from 15-25% for different smoothing methods ( Figure 5 and Table A1), in which 15% (averaged fitness is 103.50 for all smoothing methods) and 20% (average fitness is 103.86 for all smoothing methods) are the two thresholds with the lowest fitness. For all smoothing methods, the fitness gradually increases as the DT extraction threshold increases and decreases relative to their corresponding best thresholds, showing a typical 'U' shape ( Figure 5). However, the fitness variations among different thresholds for O-WS are the smallest (SD = 17.25), and the smaller fitness variation compared with other smoothing methods correspond well with the smaller R 2 and RMSE variation. As for the Zhang method for SOS extraction, the fitness for all smoothing methods are significantly lower than DTs, and the average fitness for Zhang (64.43) is 40.07 lower than that for the best threshold of DT ( Figure 5 and Table A1), with both reduced RMSE and increased R 2 ( Figure 6). Moreover, the four smoothing methods with the optimal parameters all achieve a fitness below 50

Inter-Comparison of SOS Extraction Using Different Smoothing Methods
The default parameters of all the smoothing methods result in an overestimation of SOS regarding DT in station 1, 2 and 7, which is especially severe for station 1 and 2, with an overestimation of SOS of~20 days (Figure 2b). However, for Zhang method, almost all the smoothing methods with the default parameter result in prediction of SOS close to the observation level except for WS, which shows a pervasive underestimation of SOS of 35.79 ± 16.51 for all the stations and years (Figure 2c). It is interesting to note that there are four points (marked with red circles in Figure 2b,c) with large underestimation (up to~80 days) of SOS for both DT and Zhang SOS extraction methods, and most of the smoothing methods show this underestimation with varying degrees, in which AG and DL have the greatest underestimations. As for the SOS extraction accuracy, most of the smoothing methods achieve a RMSE less than 16 days for both DT and Zhang methods (Tables 2 and 3), which is acceptable according to the temporal interval of the NDVI data [21], but WS with default parameters results in a RMSE that is even larger than 2 times the temporal interval of the NDVI data (Table 3). It is worth noting that even with the default parameters, all the smoothing methods have different performances for different SOS extraction methods, and for most of the smoothing methods, the fitness for Zhang is significantly lower than that for DT, with an average fitness 75.36 lower for Zhang compared with that for DT. For the two phenology extraction methods, the smoothing methods that lead to the highest (HANTS) and lowest (WS) accuracy are the same, in which WS has much greater RMSE for Zhang method, while the R 2 of HANTS SOS almost doubled for Zhang method compared with DT.
The optimal parameters improve the SOS extraction accuracies for all the smoothing methods compared with the default parameters, and the fitness decrease from 84. 38 (Tables 2 and 3), and the average fitness of the four optimal smoothing methods for DT and Zhang are 89.04 (decreased by 43.24% compared with the default parameters) and 41.68 (decreased by 77.49% compared with the default parameters), respectively. WS has the greatest improvement in SOS extraction accuracy for both DT and Zhang after adopting the optimal parameters, with a RMSE decrease and R 2 increase of 27.81 (3.82) and 24% (19%) for Zhang (DT), respectively. Moreover, the underestimation of SOS of WS for Zhang method has been significantly improved after using the optimal parameters (Figures 2c and 4b).
For the DT methods with nine thresholds for SOS extraction, O-CW, AG and DL result in similar fitness that are above 150, while the fitness for O-SG are generally lower than the 3 methods, but the values are still all over 100 ( Figure 5). On the contrary, O-WS and O-HANTS has much smaller fitness compared with the other 4 smoothing methods for the 9 DTs, with an average fitness value less than 90 for the 9 DTs (Table A1). The low fitness of O-HANTS and O-WS can be attributed to both higher R 2 and smaller RMSE compared with other smoothing methods (average R 2 and RMSE are 0.23 (0.23) and 13.75 (19.30) for O-WS (O-HANTS)). However, O-HANTS has higher R 2 compared with O-WS for the best threshold, while O-WS shows small RMSE (RMSE < 12) with a wide range of thresholds (15%-40%), leading to a close fitness for these thresholds. The best thresholds of DT for grassland SOS ranges from 15-25% for different smoothing methods ( Figure 5 and Table A1), in which 15% (averaged fitness is 103.50 for all smoothing methods) and 20% (average fitness is 103.86 for all smoothing methods) are the two thresholds with the lowest fitness. For all smoothing methods, the fitness gradually increases as the DT extraction threshold increases and decreases relative to their corresponding best thresholds, showing a typical 'U' shape ( Figure 5). However, the fitness variations among different thresholds for O-WS are the smallest (SD = 17.25), and the smaller fitness variation compared with other smoothing methods correspond well with the smaller R 2 and RMSE variation. As for the Zhang method for SOS extraction, the fitness for all smoothing methods are significantly lower than DTs, and the average fitness for Zhang (64.43) is 40.07 lower than that for the best threshold of DT ( Figure 5 and Table A1), with both reduced RMSE and increased R 2 ( Figure 6). Moreover, the four smoothing methods with the optimal parameters all achieve a fitness below 50 (R 2 > 0.25 and RMSE < 12.5) ( Figure 6 and Table A1). Compared with the best DT, the SOS extracted by Zhang are closer and more consistent with the observed SOS and less overestimation or underestimation can be found (Figure 6), illustrating a more stable and accurate SOS extracted by Zhang compared with DT. O-HANTS and O-WS are the two methods with the lowest fitness for Zhang, and AG and DL are the two methods with the highest fitness. Taking all the SOS extraction methods (9 DTs and Zhang method) into account, O-HANTS has the lowest fitness of 33.52, followed by O-WS, which is 4.93 higher than that of O-HANTS (Table A1).
Remote Sens. 2018, 5, x FOR PEER REVIEW 11 of 26 (R 2 > 0.25 and RMSE < 12.5) ( Figure 6 and Table A1). Compared with the best DT, the SOS extracted by Zhang are closer and more consistent with the observed SOS and less overestimation or underestimation can be found (Figure 6), illustrating a more stable and accurate SOS extracted by Zhang compared with DT. O-HANTS and O-WS are the two methods with the lowest fitness for Zhang, and AG and DL are the two methods with the highest fitness. Taking all the SOS extraction methods (9 DTs and Zhang method) into account, O-HANTS has the lowest fitness of 33.52, followed by O-WS, which is 4.93 higher than that of O-HANTS (Table A1). Figure 7 shows the difference between the SOS results of the different smoothing methods and the extraction methods. Statistically significant differences are found for all the smoothing methods except for AG and DL, and SG and HANTS. However, there are statistically significant differences for different extraction methods. Interestingly, the results of F-test of analysis of variance show that the interaction between the smoothing method and the extraction method is extremely significant (p = 0.000).      Figure 7 shows the difference between the SOS results of the different smoothing methods and the extraction methods. Statistically significant differences are found for all the smoothing methods except for AG and DL, and SG and HANTS. However, there are statistically significant differences for different extraction methods. Interestingly, the results of F-test of analysis of variance show that the interaction between the smoothing method and the extraction method is extremely significant (p = 0.000).

Parameter Sensitivity of Different Smoothing Methods
The fitness achieved by different smoothing parameters of different smoothing methods for DT and Zhang SOS extraction methods are presented in Figures 8 and 9. We only present the DT results with the best thresholds, and for each smoothing methods the threshold achieving the lowest fitness is used. On general, the changing trends of fitness for different parameters of the smoothing methods are very similar for DT and Zhang methods (Figures 8 and 9), but Zhang has lower fitness with 94.81% of the parameter values compared with DT, resulting in the average fitness of Zhang is 5609.66 lower than DT. As for DT, the fitness of CW and SG with different parameters are all over 100 (equal to the situation with RMSE = 16 and R 2 = 0.16 (R = 0.4)), and although HANTS and WS can achieve fitness below 100 or even below 50 with the right parameterization, the proportion of such parameters are less than 30% for HANTS and 8% for WS of all parameters (Figure 8). On the contrary, most of the smoothing methods (except for WS) can result in a fitness less than 100 for Zhang with most of the parameter values, especially for CW, the fitness for all parameter values for which are less than 60 ( Figure 9). No matter for DT or Zhang methods, the fitness variation of CW and SG with different parameters are much smaller than HANTS and WS, where the fitness SD for CW is 6.15 and 1.38 for DT and Zhang, respectively, which is the smallest among all methods. Oppositely, HANTS has the largest fitness range, with a fitness SD of 3142.98 for DT and 6697.88 for Zhang method, followed by WS with a fitness SD of 3126.62 and 208.71 for DT and Zhang method, respectively.
Although the parameters have different effects for different phenology extraction methods (Figures 8-10), they show some regular patterns in common. As for CW, different parameter values show small fitness variation for both DT and Zhang, and a lower fitness can be achieve with larger r and smaller fet (Figure 8), which also correspond well with the fact that r and fet are the two dominant factors for SOS extraction accuracy, and r has the most dominant impact for most cases (Figure 9). The two parameters (m and d) of SG have roughly the same contribution (45.48% from m vs. 54.52% from d on average) for the fitness value of SOS extraction (Figure 10b), and a lower fitness appears in the regions with lower m and larger d (Figures 8 and 9). WS and HANTS can have lower fitness compared with CW and SG, but only within a narrow parameter range (Figures 8 and 9). Furthermore, the parameter ranges of HANTS and WS that can result in high SOS extraction accuracy are very close for DT and Zhang method. Regarding WS, a low fitness can be achieved with parameter d set to 3 for DT or parameter d set to 3 or 4 for Zhang method, and generally as parameter λ decreases

Parameter Sensitivity of Different Smoothing Methods
The fitness achieved by different smoothing parameters of different smoothing methods for DT and Zhang SOS extraction methods are presented in Figures 8 and 9. We only present the DT results with the best thresholds, and for each smoothing methods the threshold achieving the lowest fitness is used. On general, the changing trends of fitness for different parameters of the smoothing methods are very similar for DT and Zhang methods (Figures 8 and 9), but Zhang has lower fitness with 94.81% of the parameter values compared with DT, resulting in the average fitness of Zhang is 5609.66 lower than DT. As for DT, the fitness of CW and SG with different parameters are all over 100 (equal to the situation with RMSE = 16 and R 2 = 0.16 (R = 0.4)), and although HANTS and WS can achieve fitness below 100 or even below 50 with the right parameterization, the proportion of such parameters are less than 30% for HANTS and 8% for WS of all parameters (Figure 8). On the contrary, most of the smoothing methods (except for WS) can result in a fitness less than 100 for Zhang with most of the parameter values, especially for CW, the fitness for all parameter values for which are less than 60 ( Figure 9). No matter for DT or Zhang methods, the fitness variation of CW and SG with different parameters are much smaller than HANTS and WS, where the fitness SD for CW is 6.15 and 1.38 for DT and Zhang, respectively, which is the smallest among all methods. Oppositely, HANTS has the largest fitness range, with a fitness SD of 3142.98 for DT and 6697.88 for Zhang method, followed by WS with a fitness SD of 3126.62 and 208.71 for DT and Zhang method, respectively.
Although the parameters have different effects for different phenology extraction methods (Figures 8-10), they show some regular patterns in common. As for CW, different parameter values show small fitness variation for both DT and Zhang, and a lower fitness can be achieve with larger r and smaller fet (Figure 8), which also correspond well with the fact that r and fet are the two dominant factors for SOS extraction accuracy, and r has the most dominant impact for most cases (Figure 9). The two parameters (m and d) of SG have roughly the same contribution (45.48% from m vs. 54.52% from d on average) for the fitness value of SOS extraction (Figure 10b), and a lower fitness appears in the regions with lower m and larger d (Figures 8 and 9). WS and HANTS can have lower fitness compared with CW and SG, but only within a narrow parameter range (Figures 8 and 9). Furthermore, the parameter ranges of HANTS and WS that can result in high SOS extraction accuracy are very close for DT and Zhang method. Regarding WS, a low fitness can be achieved with parameter d set to 3 for DT or parameter d set to 3 or 4 for Zhang method, and generally as parameter λ decreases the fitness decrease. It is surprising to find that the WS parameter contribution to the fitness varies much for different thresholds of DT, and as the threshold increases the contribution of parameter d increases (from 15.75 to 94.81) (Figure 10c). In terms of HANTS, a nf set to 3 or 4 can lead to lower fitness, and a lower fet can bring lower fitness (Figures 8 and 9). It should be noticed that when nf is set to 2, the fitness increases dramatically for both DT and Zhang (increases by 49 times for DT and 529 times for Zhang method), suggesting that nf = 2 is a bad option for vegetation growth smoothing.  (Figure 10c). In terms of HANTS, a nf set to 3 or 4 can lead to lower fitness, and a lower fet can bring lower fitness (Figures 8 and 9). It should be noticed that when nf is set to 2, the fitness increases dramatically for both DT and Zhang (increases by 49 times for DT and 529 times for Zhang method), suggesting that nf = 2 is a bad option for vegetation growth smoothing.

The Impact of the Choice of Smoothing Methods on the Accuracy of SOS Extraction
In this paper, it was found that for a certain SOS extraction method, different choices of smoothing methods can lead to large differences in phenology extraction accuracy, which is consistent with the conclusion of Atkinson et al. [32] that is based on four smoothing methods for the phenology extraction of the major vegetation types in India. This phenomenon also emphasizes that different combinations of phenology extraction methods and smoothing methods will lead to different accuracies of phenology extraction, but previous phenology extraction researches applied the ensemble mean of several methods as the main conclusion, and the individual method as the

The Impact of the Choice of Smoothing Methods on the Accuracy of SOS Extraction
In this paper, it was found that for a certain SOS extraction method, different choices of smoothing methods can lead to large differences in phenology extraction accuracy, which is consistent with the conclusion of Atkinson et al. [32] that is based on four smoothing methods for the phenology extraction of the major vegetation types in India. This phenomenon also emphasizes that different combinations of phenology extraction methods and smoothing methods will lead to different accuracies of phenology extraction, but previous phenology extraction researches applied the ensemble mean of several methods as the main conclusion, and the individual method as the complementary analysis [61,62], which barely consider the impact of the interaction of phenology extraction methods and smoothing methods. We found that the values of smoothing parameters have a great impact on the phenology extraction results, which is similar to the results of Cai et al. [63] who found that smoothing methods can have different curve smoothing performances and curve shapes with different parameters, which can be explained that the valuing of the smoothing parameters affects the curve smoothness as well as the value ranges of the smoothed curves (Figures 2 and 4).
Our study showed that for most of the smoothing methods, the default parameters can bring acceptable SOS extraction results, while the default parameters for WS showed poor SOS extraction performance regarding both DT and Zhang methods, which may be because the default parameters of WS in Atkinson et al. [32] are set oriented for multiple vegetation types rather than for grassland. Besides, our results showed that WS can have high SOS extraction accuracy for a large range (10% to 45%) of the threshold of DT method ( Figure 5) with specific parameter optimizations, although a high SOS extraction accuracy is preferable, a high accuracy achieved at an unsuitable threshold suggests that the shape of WS may be unstable and can change dramatically with different parameters, which can also be confirmed that the contributions of the WS parameters vary across different DT thresholds (Figure 10c). The findings in our study suggest that when using WS to extract vegetation phenology, careful attention should be paid to its parameters, and a parameter optimization may be required to achieve a satisfactory result. The two non-parameter smoothing method (AG and DL) have always resulted in lower SOS extraction accuracies than other methods, which is because AG and DL adopt Gaussian and Logistic function to fit the NDVI curve, and the fitted curves are susceptible to low-value noise, resulting in changing curve shape (Figure 11b-d) and leading to a lower SOS extraction fitness. HANTS was found to have the highest SOS extraction accuracy no matter with the optimal parameters or the default parameters, which is because HANTS can maintain the annual minimum (Figures 2a  and 4a), and the slope of the curve around SOS is stable (Figure 3). The default parameter set and the optimal parameter sets for HANTS are very close (Table 2), suggest the default parameter set for HANTS is a good choice for grassland SOS extraction.
Previous studies evaluate the performances of smoothing methods mostly aiming at their denoising effect against the simulated noise [54,64]. However, by comparing different smoothing methods and their SOS extraction accuracy, we found that the performance of curve denoising does not definitely match the accuracy of SOS extraction. For example, the errors introduced by CW and SG during their preprocessing of the continuous noise destroy the correct shape of the curve and cause new noise in some cases, reducing the accuracy of DT-based SOS extraction (Figure 3c). On the other hand, the smoothing result of the O-WS introduces many extreme low values in the right half curve using the Zhang method ( Figure A2), but the growth rate of the grassland is maintained well in the first half of the smoothed curves, leading to a high SOS extraction accuracy ( Figure A2). Overall, based on our findings, the recommended smoothing method for grassland SOS extraction is HANTS, not only for its highest SOS extraction accuracies with both default and optimal parameters, but its stable performance for different SOS extraction methods as well. Although HANTS has the large variation of SOS extraction accuracy with different parameters, it is mainly due to the magnificent fitness values when parameter nf is set to 2, and after excluding value 2 from the nf value range, the fitness variation of HANTS (SD = 64.00) drops dramatically to the CW and SG level.

The Impact of Smoothing Method Parameters on the Accuracy of SOS Extraction
Different smoothing methods show different parameter sensitivities, and CW was found to have a lesser impact on parameters than the other methods. It can be explained that parameter r is the dominant parameter impacting the SOS extraction accuracy of CW (Figure 10a), and as the searching radius of the local maximum and minimum, parameter r is the important parameter of preserving original maximum and minimum. However, the variation of the minima in left half of the curve with different r is small for single growing season vegetation (e.g., grassland). Furthermore, the three-point weight strategy that CW adopts can preserve the whole shape of the original curve with different local minimum and maximum [34]. Therefore, no matter for DT and Zhang, the SD of O-CW with different parameters is far smaller than others. As for SG, our results are consistent with Chen et al.'s [35] findings that smaller value m and a larger d leading to a better smoothing time-series. Atkinson et al. [32] concluded that the use of a smaller parameter λ of WS will produce larger errors, a λ = 15 is suitable for vegetation phenology extraction. However, in this paper, the SOS extraction accuracy with a λ = 15 is much lower than with a smaller λ (e.g., λ = 2). It may be because that Atkinson et al. [32] carried out their study based on a variety of vegetation, and since the growth curves of different vegetation differ greatly [65], a smoother curve can to blur the differences among vegetation to obtain an overall high accuracy of SOS extraction. For the SOS extraction of grassland in our study, more accurate characterization may be required, and when the smoothed curves are too smooth, the under-fitted results may lead to errors. There are no objective rules for determining parameters for HANTS [36], and the nf is believed as a key parameter for HANTS [54]. Based on our findings, parameter nf of HANTS was recommended to set between 3-4, which is close to the optimizations of previous studies [36,54,66]. As parameter nf represents the number of frequencies of the Fourier components in HANTS [36], a large value nf makes HANTS easily affected by noise and prone to have multi-peaks in the smoothed curves. Besides, we found a large nf can result in a failure of smoothing under some of the condition with a large number of noise during the simulations, so a smaller nf is more preferable for vegetation indices time-series smoothing, especially at large scales.
Here, suggestions about the parameter optimizations of different smoothing methods for grassland SOS extraction are drawn based on our findings. As for CW, parameter setting has little impact on the SOS extraction accuracy, the default parameter setting is good enough to achieve an acceptable or good SOS extraction result. However, if a parameter optimization is required, parameter r is the key parameter that need to be taken into consideration first, followed by parameter fet. Generally, a smaller r and fet is recommended for grassland SOS extraction; the two parameters of SG have the same influence on SOS extraction, where parameter d is not suggested to be set to a too small value; for WS method, parameter d is recommended to set between 3 and 4, and a smaller λ can bring higher SOS extraction accuracy; Similar to WS, parameter nf of HANTS is suggested to set between 3 and 4 for a higher SOS extraction accuracy, and a smaller fet is also suggested.

The Impacting Factors on SOS Extraction Accuracies for Phenology Extraction Methods
In this study, the best threshold of DT for grassland SOS extraction was found to be between 15-25%, which is in line with the widely adopted 20% threshold in previous studies [21,26]. However, the best thresholds vary cross different smoothing methods, and even with the best threshold, the SOS extraction accuracy of DT is commonly lower than Zhang method. According to Sections 4.1 and 4.2, we attribute the varied thresholds and low accuracy of DT to the fact that the performance of DT is greatly affected by the fluctuations of the minimum value in the NDVI curves (Figures 2a and 4a). To solve this problem, Yu et al. [21] manually fixed the minimum value of NDVI, and the DT with a fixed minimum value (hereinafter referred to as DT-fix) was used for grassland SOS extraction. To further validate the effects of low-value noise on DT SOS extraction, we followed the concept of Yu et al. [21] to fix the NDVI minimum value, and to simplify the simulating process, the minimum NDVI value was fixed to 0.05 to represent a bare soil condition according to Hird and McDermid [67], and the result from DT-fix were compared with that from DT and Zhang methods. It is easy to find that DT-fix substantially improves SOS extraction accuracy from DT, with both reduced RMSE (decreased by 3.13 for default parameter and 2.44 for optimal parameter on average) and increased R 2 (increased by 12.97 for default parameter and 11.79 for optimal parameter on average) (Tables 2 and 4), resulting in 49.72% and 23.70% fitness reduction for default parameters and optimal parameters, respectively. Moreover, for the most cases, DT-fix has even achieved a lower fitness (higher SOS extraction accuracy) compared with Zhang method for both default parameters and optimal parameters, with a fitness 71.65% and 7.54% lower on average (Tables 3-5). The best thresholds of DT-fix for different smoothing methods are all 35%, and the higher threshold compared with 20% may be due to that 0.05 is much lower than annual NDVI minimum of grassland on QTP. However, the uniform threshold highlighted that DT-fix indeed suffer less from the adverse effects from the fluctuations of NDVI minima. The DT-fix results showed that with the fixed NDVI minimum, the negative effects from NDVI minima fluctuations are mitigated to a great extent, and DT-fix has higher accuracy for SOS extraction compared with DT and Zhang methods. In our study, the Zhang method showed superior SOS extraction accuracies for grasslands compared with DT, which is consistent with a previous study that based on winter wheat [57]. This could be explained by the idea that the Zhang method is to fit the curve with the DL function and calculate the local maximum value of the curvature transformation rate of the fitted curve [15], so the shape of the first half of the curve, especially the shape of the first half curve around SOS, is the main factor affecting the extraction accuracy of SOS, while the fluctuations of NDVI minimum values have less impact on the Zhang method compared with DT [68], and the better performance of Zhang may be due to that the minimum of first half of the curve is more susceptible by noise than shape of the first half curve around SOS after smoothing in the most instances. This can also be confirmed by the curves of the four red-circled points in Figure 2c, and the NDVI curves smoothed by different smoothing methods of the 4 points are shown in Figure 11. In point 1, the curve shapes of CW and SG deviates farther from those of other methods (Figure 11a), while at the remaining three points (Figure 11b-d), similar deviations can be found for AG and DL, and these deviations of the curve shapes correspond well with the differences in the SOS extraction accuracies.

Applicability of Different Smoothing Methods
In this study, the grassland SOS extracted using different smoothing methods and different phenology extraction methods were compared, and the results can provide some guiding insights for the choice of smoothing methods and phenology extraction methods under different scenarios. Firstly, when there are no available observed data to optimize smoothing parameters, HANTS is the most recommended method due to its high SOS extraction accuracies for both DT and Zhang with the default parameters. Besides, CW and SG are also recommended because of their acceptable (for DT) and good (for Zhang) accuracies of SOS extraction and stable performance (less sensitive to parameter settings). AG, DL and WS with the default parameters is not suggested for grassland SOS extraction, mainly due to their low accuracies of SOS extraction. Secondly, if there are abundant phenology observations to optimize the smoothing parameters, HANTS and WS are more recommended due to their excellent for both DT and Zhang performance after adopting the optimal parameters.
As for the phenology extraction methods, the Zhang method has obvious advantages over DT, and is more suggested for grassland SOS extraction. However, DT can have a robust and accurate performance of SOS extraction by fixing the minimum value, with a SOS extraction accuracy even higher than the Zhang method. The threshold of DT-fix in suggestion for grassland SOS extraction is 35% based on our results.

Conclusions
In this paper, based on MOD13Q1 NDVI time-series and observed SOS records of 8 agrometeorological stations, we compared six smoothing methods for grassland SOS extraction on the Qinghai-Tibetan Plateau during 2001-2011. We found that the bad denoising performance is not in line with the low SOS extraction accuracy (e.g., O-WS achieved a low fitness for Zhang method even with extremely low values in the smoothed curve). Different smoothing methods show different parameter sensitivities, and the optimal parameters can improve the accuracy of SOS extraction. In addition, the optimal parameters are different for different extraction methods, the average fitness of the four optimal smoothing methods for DT and Zhang are decreased by 43.24% and 77.49% compared with the default parameters, respectively. For the 6 smoothing methods, HANTS has lowest fitness for Zhang (fitness = 33.52 with parameter optimization and fitness = 35.05 without parameter optimization) and the denoising ability of HANTS are all better than other methods. The phenology extraction method has a greater impact on accuracy than the smoothing method and the main influencing factor of DT and Zhang are the stability of the annual minimum and curve shape near SOS, respectively. Zhang is better and more stable than DT method for all smoothing methods. However, after setting the minimum value of the NDVI curve to eliminate the fluctuation error of the annual NDVI minimum value, with a suitable threshold, the DT-fix method has less difference in fitness for the default parameters or optimal parameters for all smoothing methods and can achieve better SOS extraction results than Zhang, with a fitness decreased by 71.65% for default parameters and 7.54% for optimal parameters on average.

Acknowledgments:
The authors would like to acknowledge China Meteorological Administration for providing ground observed phenology data.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.