## 1 Introduction  
Accurate time-series of vegetation indices from satellite remote sensing are vital for long term vegetation monitoring, especially for studies of vegetation phenology [1] and carbon balance [2; 3]. However, the signals received by the satellite sensors are affected by noise due to geometric misregistration, anisotropic reflectance effects, electronic errors, artifacts due to data resampling, atmosphere, and clouds [4]. In order to reconstruct the vegetation seasonal growth course from noisy satellite signals a multitude of time-series processing methods have been used. Some early methods are e.g., the best index slope extraction [5], the repeated moving-window median filter [6; 7; 8], harmonic functions [9] and median filters [10]. More recently and widely used methods are e.g., Savitzky-Golay filtering [11], least-squares fits to asymmetric Gaussian functions [12] and double logistic functions [13; 14; 15; 16], and variations of spline smoothing [17; 18] including the Whittaker smoother [19], and wavelet transforms [20; 21; 22]. The processing methods can be divided into two broad categories: (1) methods that enable precise capturing of short-term variations during the growing season, generally by employing different types of statistical filters [5; 6; 10]; and (2), methods that fit mathematical functions (sinusoidal, logistic, Gaussian etc. [12; 13; 14; 15; 16; 20; 21; 22]) to part of a season, the full season, or even a sequence of seasons. We can in a loose sense label these as \"local\" and \"global\", respectively, denoting their flexibility to match to details in the time-series or to generalize across the full growing season. Spline methods [17; 18; 19] have the ability to take on both local and global character depending on a smoothing parameter. Due to the difficulty of existing physiological process-based phenology models, e.g., various budburst models using chilling or forcing temperature requirements, in accurately simulating phenology from climate data [23], time-series of smoothed vegetation indices measured by satellite are an important data resource for estimating phenology [24; 25]. Many studies have employed remote sensing time-series to extract and monitor the timing of plant phenological events [13; 15; 26; 27; 28; 29]. Since these studies have used different smoothing methods to estimate phenology, it is hard to draw any firm conclusions about which methods are the most suitable. [PERSON], et al. [30] compared 10 different phenology retrieval methods applied to an Advanced Very High Resolution Radiometer (AVHRR) NDVI dataset in North America and concluded that although smoothing methods applied to satellite NDVI could generally capture vegetation phenology, no method was universally better than the others. [PERSON], et al. [31] and [PERSON], et al. [32] concluded that the performance of smoothing methods in estimating vegetation phenology varies spatially and temporally, and that, because of bias and random errors due to cloud influence, no single method exhibited a consistently superior performance. [PERSON], et al. [33] argued that the reason for not finding a single best method could be different definitions of phenology parameters. [PERSON] and [PERSON] [34] highlighted that the performance of different smoothing methods depends on both land surface condition and the clear sky identification approach adopted. Many studies have investigated the performance of smoothing methods in reconstructing consecutive time-series of vegetation indices, leaf area index and fractional absorbed photosynthetically active radiation (fAPAR) [31; 35; 36; 37; 38; 39]. Besides vegetation growth time-series from satellite remote sensing, ground-based optical measurements, e.g., Phenocam [40] and multispectral sensors [41; 42], have been successfully applied in vegetation phenology studies for different biomes [29] and for linking to satellite remote sensing measurements [43; 44].