A New Equation for Deriving Vegetation Phenophase from Time Series of Leaf Area Index ( LAI ) Data

Accurately modeling the land surface phenology based on satellite data is very important to the study of vegetation ecological dynamics and the related ecosystem process. In this study, we developed a Sigmoid curve (S-curve) function by integrating an asymmetric Gaussian function and a logistic function to fit the leaf area index (LAI) curve. We applied the resulting asymptotic lines and the curvature extrema to derive the vegetation phenophases of germination, green-up, maturity, senescence, defoliation and dormancy. The new proposed S-curve function has been tested in a specific area (Shangdong Province, China), characterized by a specific pattern in leaf area index (LAI) time course due to the dominant presence of crops. The function has not yet received any global testing. The identified phenophases were validated against measurement stations in Shandong Province. (i) From the site-scale comparison, we find that the detected phenophases using the S-curve (SC) algorithm are more consistent with the observations than using the logistic (LC) algorithm and the asymmetric Gaussian (AG) algorithm, especially for the germination and dormancy. The phenological recognition rates (PRRs) of the SC algorithm are obviously higher than those of two other algorithms. The S-curve function fits the LAI curve much better than the logistic function and asymmetric Gaussian function; (ii) The retrieval results OPEN ACCESS Remote Sens. 2014, 6 5651 of the SC algorithm are reliable and in close proximity to the green-up observed data whether using the AVHRR LAI or the improved MODIS LAI. Three inversion algorithms shows the retrieval results based on AVHRR LAI are all later than based on improved MODIS LAI. The bias statistics reveal that the retrieval results based on the AVHRR LAI datasets are more reasonable than based on the improved MODIS LAI datasets. Overall, the S-curve algorithm has the advantage of deriving vegetation phenophases across time and space as compared to the LC algorithm and the AG algorithm. With the SC algorithm, the vegetation phenophases can be extracted more effectively.


Introduction
Vegetation phenology is the seasonal phenomenon of vegetation, which is influenced by the environment, especially long-term temperature changes [1].Phenological events include germination, foliage expansion, foliage maturation, fall foliage coloration, defoliation and dormancy.Phenology not only reflects the seasonal alternation, but also the adaptability of vegetation to environmental conditions.With the changing of the global climate over nearly half a century, the phenology of vegetation has adjusted to ensure survival and reproduction [2], and these phenological changes are the most sensitive indicator of climate change [3][4][5][6].Thus, understanding changes in vegetation phenology is very important to the study of vegetation ecological dynamics, comprehending the influence of climate change on vegetation, forecasting the trend of climate change and identifying ways to slow climate change [7].
So far, three main types of approaches can help people identify the vegetation phenophases [8,9], which are phenological observations on the ground for individual plants or tree stands [4,10], derivative phenophases based on remote senesing data [2,6,8,[11][12][13][14][15][16][17], and plant phenology models based on climate factors (e.g., temperature, light and soil moisture, etc.) [18][19][20].Among these approaches, the derivative phenology using time series of remote sensing data, which have continuous coverage of time and space, has been approved as a promising and reliable measure at the large spatial scales [8,21], and widely applied in vegetation phenology studies.These applications are mainly based on three datasets: MODIS (MODerate-resolution Imaging Spectroradiometer), AVHRR (Advanced Very High Resolution Radiometer) and SPOT-VGT (Satellite Pour l'Observation de la Terre-Vegetation).The variables used to detect the vegetation phenophases mainly include the vegetation index (VI) and leaf area index (LAI).Thus, several vegetation phenology retrieval methodologies have been designed and include the thresholds, the moving average, and the fitting function of these variables.
The threshold methods contain static and dynamic approaches, and are often applied after filtering the VI or LAI data to reduce noises and produce a smoother time series.The filter approaches involve the harmonic analysis of time series (HANTS) algorithms [5,22,23], the Savitzky-Golay filtered algorithm [6,16], the Fourier filter [21,24], the Whittaker smoother [11] and the six-time polynomials [5,25,26].The static threshold method identifies the start or end of the vegetation growing season by setting a constant VI or LAI value [27].The static threshold method is limited in that it does not always apply to the entire region.In contrast, the dynamic threshold extracts the vegetation phenology by using the varying VI or LAI value obtained from various threshold algorithms [11,15,26,[28][29][30].Although the advantage of evaluating vegetation phenology using the dynamic threshold is obvious, the accuracy of the retrieved results strongly depends on the threshold algorithms.
The moving average method identifies the start or end of the vegetation growing season based on the date that corresponds to the intersection between the time series of the VI or LAI up or down curve and the moving average curve [31,32].The accuracy of the retrieved results and the inefficacy in recognizing phenological changes are closely related to the moving window size.
The fitting function method refers to the ecological function and the corresponding phenophases extraction method.This retrieval method first applies the ecological function to fit the VI or LAI data to get a fitted curve, then uses the phenophases extraction method to detect the vegetation phenology based on the fitted curve.The ecological functions used frequently include the asymmetric Gaussian function [21,33], the double logistic function [21,34] and the piecewise logistic function [13,[35][36][37][38][39].The phenophases extraction methods usually contain the change rate or the curvature of the fitted curve.The fitting function method describes the trends of the variation of vegetation phenology in time and space more effectively than the threshold and the moving average method, and it does not require thresholds or empirical constraint conditions.Therefore, the fitting function method is broadly applicable to deriving vegetation phenology; however, its ability to recognize phenology is lower than that of the other methods due to the ecological function requiring the VI or LAI curve possesses the Sigmoid feature.
To date, these phenology retrieval methods have focused on detecting the start and end of the vegetation growing season [5,11,15,[26][27][28][29][30][31].The associated retrieval results are inconsistent [21].Based on the ecological significance of the algorithms, the retrieval results of the fitting function method are more reasonable.Thus, the piecewise logistic function method has been applied to the MODIS global land cover dynamics product (MLCD) [36].However, the failure rate of the MLCD retrieval results is also high [40].Therefore, the logistic function method must be further improved to increase the accuracy of the retrieval results.With regard to the asymmetric Gaussian and the double logistic methods, they did not perform well for areas with more than one growing season per year [21].Considering the conditions described above, we developed a more universal phenology inversion algorithm based on the asymmetric Gaussian function and the logistic function.With this algorithm, the vegetation phenophases can be extracted more effectively.

Site Descriptions
In the study, we selected 10 national agriculture and meteorology observation stations in China.They are the Heze, Liaocheng, Jining, Taian, Huimin, Zibo, Linyi, Weifang, Laiyang and Wendeng sites (Figure 1a).These stations not only measure the agriculture and meteorology variables, but also provide the information on vegetation phenology.We thus call it a "phenology observation station" here.The site descriptions are shown in Table 1.

Figure 1.
The terrain in Shandong Province (a), the spatial distribution of the phenology stations (a) and the spatial distribution of the vegetation types (b).In Figure 1b, the abbreviation NET is temperate needle-leaf evergreen trees; BDT is temperate broad-leaf deciduous trees; BDS is temperate broad-leaf deciduous shrubs.
The selected phenology stations are all situated in Shandong Province, which is located on the eastern coast of China and in the downstream area of the Yellow River (114.32-122.72°E,34.37-38.38°N).The central and southern regions of the province are mountainous and hilly.The eastern peninsula of Shandong Province is a hilly area, and the western and northern areas are within the Yellow River alluvial plain (Figure 1a).The climate is moderate and is characterized as warm temperate continental monsoonal.The seasons are easily distinguished here, and the rainy season is primarily in the summer.The multi-year average air temperature is 11-14 °C, and the annual average precipitation is approximately 500-1100 mm.The vegetation type across the province is mainly crops.The mountainous and hilly areas are dominantly covered by grass, temperate needle-leaf forest and broadleaf forest (Figure 1b).

LAI Datasets
The LAI datasets used in this paper contain the improved MODIS LAI and AVHRR LAI-3g (uses of the term "AVHRR LAI" occurring in the following text all refer to "AVHRR LAI-3g", unless specified otherwise).The former is an improved version of the standard MODIS LAI product (MOD15A2), which provides 8-day global LAI data from 2000 to 2009 at a 1 km spatial resolution [41].The latter was generated from the AVHRR GIMMS (global inventory modeling and mapping studies) NDVI-3g (normalized difference vegetation index, third generation) data by using the Feed-Forward Neural Network algorithm based on the training data of improved MODIS LAI [42].The AVHRR LAI dataset has a spatial resolution of 1/12 degree (nearly 8 km) and a 15-day interval, and spans the period from 1982 to 2009.The GIMMS NDVI-3g dataset improved AVHRR data quality in the northern parts of the world by using improved calibration procedures and calibrated using Sea-viewing Wide Field-of-view Sensor (SeaWiFS) data [8], and has been corrected to minimize various deleterious effects, such as calibration loss, orbital drift and volcanic eruptions, and has been verified using stable desert control points [43,44].The improved MODIS LAI dataset greatly removed the noises due to the presence of cloud, seasonal snow cover, the instrument problems, and the uncertainties of retrieval algorithm [41].Based on the reliable GIMMS NDVI and improved MODIS LAI, the quality of derived AVHRR LAI is good and reliable.However, the two LAI datasets still contain some undetected noises from cloud contamination, atmospheric variability and bi-directional effects.To further eliminate these noises, both two LAI datasets are first subjected to the iterative Savitzky-Golay filtered algorithm.We treated the filtered improved MODIS LAI dataset and the AVHRR LAI dataset as the "ground truth," although some uncertainties still exist in these two LAI datasets, which are limited by other issues, e.g., uncertainties in land cover types, etc.To scale the LAI temporal frequency down from 15-day or 8-day to 1-day, we used the cubic spline interpolation, which has a superior smoothing feature and effectively reflects the signals of the original LAI data to interpolate the missing data [6,45,46].In the land cover types with well-defined phenologies, the cubic spline algorithm performed reasonably well for temporal interpolation [46].Thus the noise caused by interpolation with cubic spline algorithm can be ignored in this study.For the purpose of minimizing the background noisy of LAI data, which can be caused by deserts, sparsely vegetated areas or other factors, the pixels with an LAI value of less than 0.2 were regarded as invalid [45].Finally, based on the coordinates of phenology observation stations described in the previous section, we used the proximal sampling to collect the corresponding LAI data from the improved MODIS LAI and AVHRR LAI at each station during the study period.
Terrain information for Shandong Province (Figure 1a) was obtained from the NASA shuttle radar topographic mission (SRTM)-produced digital elevation data (DEM) with a 90 m resolution (http://www.cgiar-csi.org/).The spatial distribution of the vegetation types (Figure 1b) was obtained from the study by Ran et al. (2012) [47].

LAI Curve Fitted Model (1) Model 1: Asymmetric Gaussian function
The left half of the asymmetric Gaussian function [33] after the equivalent-transformation is described below: ( ) (( ) / ) a m t a t a  and t < a 1 , y is the LAI-fitted value and t is the time in days.The parameters p and q determine the amplitude and the base level, respectively; a 1 determines the position of the maximum or minimum with respect to the independent time variable t. a 2 and a 3 determine the width and flatness (kurtosis), respectively, of the left half of the function (in most cases, a 3 = 2).
(2) Model 2: Logistic function The logistic function [13] uses the following function form: where m 2 (t) = b 1 t + b 2 , y is the LAI-fitted value and t is the time in days.The parameters p and q determine the amplitude and the base level, respectively; b 1 and b 2 are both fitted parameters.
(3) Model 3: S-curve function From Equations ( 1) and ( 2), the uniform function is as follows: where m(t) = at 2 + bt + c; y is the LAI-fitted value; t is the time in days; q is the minimum LAI value; p + q is the maximum LAI value; and a, b and c are parameters of the curve and control the position of the inflection points as well as the slope of the curve.Two special cases are worth considering.If 3) is equivalent to the asymmetric Gaussian function of Equation ( 1).When comparing for the fit of the LAI data after the log-transformation, the quadratic polynomial (m 1 (t) and m(t)) performs better than the linear function (m 2 (t)) [45].The difference between m 1 (t) and m(t) is whether repeated roots exist in m(t).
Here, we named Equation (3) the Sigmoid curve (S-curve) function.The seasonal variation of vegetation usually has at least one peak in a given year.For example, a temperate seasonal deciduous forest has a single peak, whereas a temperate crop usually has two peaks.Thus, the use of the piecewise S-curve function to simulate the seasonal variation of LAI is necessary: ( ; , , , , ) () ( ; , , , , ) where y is the LAI fitted value; t is the time in days; y l and y r are the left side and right side LAI fitted values, respectively; T is the time in days when LAI reaches its maximum; and (p l , q l , a l , b l , c l ) and (p r , q r , a r , b r , c r ) are the sets of parameters for the left side and the right side of the S-curve function, respectively.

Vegetation Germination and Dormancy
Theoretically, an asymptotic line exists on the bottom of the S-curve (y(t) = q).Therefore, the date that corresponds to the point where the asymptotic line on the left side of the S-curve function intersects the S-curve function through a given small tolerance represents the vegetation germination (i.e., the start of the vegetation growing season) (point A in Figure 2).Similarly, the date that corresponds to the point at which the asymptotic line on the right side of the S-curve function intersects the S-curve function through a given small tolerance represents the vegetation dormancy (i.e., the end of the vegetation growing season) (point F in Figure 2).In practice, the fitted value of an S-curve in a given period usually cannot attain the q-value.Thus, the minimum of the fitted value of the S-curve serves as the corresponding asymptotic line (Figure 2).The direction convention is as follows: the direction of the left (right) side of the S-curve is toward the right (left).The point at which the asymptotic line intersects the S-curve function is calculated by the tolerance (the value in this study is set to 0.01).

Vegetation Green-Up, Maturation, Senescence and Defoliation
We use the curvature extremum methods to extract the vegetation phenophases.From Equation (3), the curvature of the S-curve function is calculated as follows: The variables y' and y'' are the first and second derivatives of the S-curve function, respectively: (2 ) (1 ) 2 (1 ) '' (1 ) where z = exp(at 2 + bt + c).
In the vegetation growth phase in spring, from the vegetation germination to the LAI maximum, the curvature k has two maxima that correspond to B and C in Figure 2. Point B is the inflection point at which the slow vegetation growth becomes rapid growth.After passing this point, the greenness index of vegetation increases quickly.Thus, this point serves as the vegetation green-up.Point C is also an inflection point that represents the most prosperous stage of vegetation growth.After passing this point, the LAI reaches the annual peak.Thus, this point is expressed as the vegetation maturation.Corresponding to the phase of vegetation growth in spring, the curvature k has two maxima during the vegetation recession in autumn (points D and E in Figure 2).Point D denotes the beginning of the vegetation senesce from the most prosperous stage.After passing this point, the LAI decreases rapidly.Thus, this point serves as the vegetation senescence or fall foliage coloration.Point E is the inflection point that marks when the vegetation senescence slows down.After passing this point, the vegetation will lose leaves; thus, this point denotes vegetation defoliation.

Extraction Process for Vegetation Phenology
Here, we defined a term "SC algorithm" which refers to the S-curve function and the phenophases extraction methodology described in previous section.Applying SC algorithm to extract phenophases follows the following steps.
Based on Figure 1b, the main vegetation type in Shandong Province is crops, followed by grass and deciduous forest (including shrubs).Because of the warm, semi-humid monsoon climate, the crops in most areas have two peaks per year.However, in a few areas, the crops have a single peak due to drought.Based on the farming work period in Shandong Province, the transition time of the double seasonal crop rotation often occurs in June: the harvest of the winter wheat often occurs in early June, and the sowing of the summer maize is concentrated in late June.Thus, the transition time of the double seasonal crop rotation was uniformly set to the end of June, which is the 181st day of the year (the black dotted line in Figure 2).In this way, the LAI data in the first and second half of the year were divided.The first dataset was used to derive the spring vegetation phenology, and the second dataset was used for the autumn phenology.For single seasonal crops and other single seasonal vegetation, the transition time was also set to the 181st day of the year.
For the LAI data in the first half of the year, the date that corresponds to the LAI maximum was calculated first.Then, the date was used to divide the LAI data into two parts: a left-side LAI dataset and a right-side LAI dataset.The left-side LAI dataset was employed to estimate the parameters (p l , q l , a l , b l , c l ) with the ordinary least square method, and detect the spring vegetation phenology (interval AC in Figure 2) with the phenophases extraction method.The same LAI data division method was applied to the LAI data in the second half of the year.However, the right-side LAI dataset was used to estimate the parameters (p r , q r , a r , b r , c r ), and identify the autumn vegetation phenology (interval DF in Figure 2).
Similar to SC algorithm, we defined the terms "AG algorithm" and "LC algorithm" as the asymmetric Gaussian function and the logistic function, respectively.Both of them employ the same phenophases extraction methodology as that in the SC algorithm.The steps for detecting phenophases with the two algorithms are similar to that of SC algorithm.
The observation sites in China's animal and plant phenology observation reports are Liaocheng, Taian and Weifang (Figure 1a).The three sites are used for validating the inversion algorithms at the single point scale.In general, each observed site is rich in its biodiversity.To reduce the difficulty of processing the phenology data of the observed sites, first major dominant and typical plant species, covering a large area at each site, were selected.Then the mean of their phenophase data was calculated for each year.We did this because we tried to reduce the deviation caused by the different spatial scales between the observed site mainly aiming at individuals of some specific plants and the pixel of the remote sensing phenology usually reflecting the compositive information of multiple vegetation types.We adopted the average strategy because we cannot calculate the weight of each dominant species.At the end, the annual phenology observation data at each site were obtained.
The sites of the national agriculture and meteorology observation stations include Heze, Liaocheng, Jining, Taian, Huimin, Zibo, Linyi, Weifang, Laiyang and Wendeng (Figure 1a).The green-up observation data, including early green-up dates and late green-up dates in spring, at these sites were obtained by digitizing the green-up points reported in the literatures [52].These green-up data were also used to validate the inversion algorithms at the single point scale.

Quantitative Evaluation of Vegetation Phenology Retrieval Methodology
The bias statistics was used to evaluate the accuracy of the phenology retrieval results based on these three inversion algorithms.It is defined as follows: where P i and O i denote the retrieval and observed values, respectively.Then, the phenology recognition rate (PRR) was used to assess the ability to recognize the vegetation phenophases using the inversion algorithms.The PRR is defined as the proportion of the number of pixels in which the inversion algorithm can recognize the phenophases compared to the number of all of the pixels involved in the calculation.Finally, to describe the efficacy of fitting the LAI data with the LAI curve fitted function, the root mean square error (RMSE) and the index of agreement (IA) were both used, which are defined below: where P i and O i denote the fitted and observed values, respectively.

Results
The most effective method for evaluating the phenology inversion algorithm is to compare the retrieval results with the observation data.We selected Liaocheng, Taian and Weifang sites in China's animal and plant phenology observation reports for the single point validation.Figure 3 shows the comparison of the phenology retrieval results with the three inversion algorithms using the AVHRR LAI datasets.The corresponding bias comparison is shown in Table 2. Figure 4 shows the comparison of the fitted LAI curves.
At the start and end of the vegetation growing season, the retrieval results of the SC algorithm are much closer to the observation data than those of using the other two algorithms (Figure 3 and Table 2).The retrieval results of the LC and AG algorithm revealed earlier dates at the start of the growing season and later dates at the end of the growing season.This difference in the retrieval results among the three inversion algorithms was shown in the fitted LAI curves (Figure 4).In the upward phase of the LAI curves (corresponding to point A in Figure 2), the S-curve function starts later than the logistic function and the asymmetric Gaussian function.However, in the downward phase of the LAI curve (corresponding to point F in Figure 2), the S-curve function begins earliest.
For the vegetation green-up and maturation in spring, the retrieval results of three inversion algorithms are all similar to the green-up observation data but dissimilar to the maturation observation data.Specifically, the green-up retrieval results of the SC algorithm are superior to those of the other two algorithms, and in close proximity to the green-up observation data.Although the results of LC algorithm were also closer to the observed, the PRRs at three sites were all the lowest (Table 3).The retrieval results of the AG algorithm all advanced lots of days in relation to the observed.The deviation between the maturation retrieval results and the observation data is obvious for all three algorithms.However, the maturation retrieval results of the AG algorithm are better than those of the other two algorithms.Based on the fitted LAI curves (Figure 4), the rapid rising phase of the S-curve function starts later than those of the other two functions.However, the slowly rising phase of the S-curve function begins earlier than that of the logistic function, and later than that of the asymmetric Gaussian function.Thus, the green-up date (point B in Figure 2) that corresponds to the curvature extremum of the S-curve function is later than those of the other two functions.Additionally, the maturation date (point C in Figure 2) that corresponds to the curvature extremum of the S-curve function is earlier than that of the logistic function, but later than that of the asymmetric Gaussian function.
Considering the vegetation senescence in autumn, the larger deviation between the retrieval results and the observed value is obvious, despite the results of AG algorithm being superior to those of the other two algorithms.For the vegetation defoliation in autumn, the accuracy of each retrieval algorithm was various at different sites.At the Taian site, for example, the results of the SC algorithm were closer to the observed than those using the other two algorithms.However, at the Weifang site, the results of the AG algorithm were the best.These differences are also reflected in the fitted LAI curves.The rapid declining phase of the S-curve function starts later than that of the logistic function, but earlier than that of the asymmetric Gaussian function (Figure 4).However, the slow declining phase of the S-curve function begins earlier than the other two functions.Thus, the senescence dates (point D in Figure 2) that correspond to the curvature extremum of the S-curve function start later than those of the logistic function, and earlier than those of the asymmetric Gaussian function.However, the defoliation dates (point E in Figure 2) that correspond to the curvature extremum of the S-curve function begin earliest.
From the comparison of the performances of the fitting LAI curves using the S-curve function, logistic function and asymmetric Gaussian function (Table 3), we can see that the RMSE of the S-curve function is the lowest among these three functions.Moreover, the IA of the S-curve function is higher than those of the other two functions.The S-curve function has a better performance in fitting the LAI curve.From Table 3, we also find the PRRs of S-curve function are the highest, suggesting that the application of the SC algorithm in the time range is superior.
In addition to the AVHRR LAI used in this study, the improved MODIS LAI was also employed.Figure 5 shows the comparison of retrieval results based on the two LAI datasets.No matter which kind of LAI dataset is used, the retrieval results using the SC algorithm were reliable and closer to the green-up observation data.The results of LC algorithm are closer to the observed values as well, just earlier than those using the SC algorithm.However, the PRRs of LC algorithm are lower than those of the other two algorithms.The accuracy of retrieved results using the AG algorithm is disappointed because its results advanced lots of days comparing with the observed values.Thus, the application of the SC algorithm is much broader than the other two algorithms across time and space.The three inversion algorithms show that the retrieval results based on AVHRR LAI were behind those based on the improved MODIS LAI.The bias statistics revealed that the retrieval results based on the AVHRR LAI dataset are more reasonable than those based on the improved MODIS LAI datasets.

Complexity of Validation of Retrieval Results with Inversion Algorithms
Conceptually, the phenophases derived from the inversion algorithms based on remote sensing data are not always in agreement with the observed phenophases.Here, six phenophases were defined according to the features of the LAI data.One phase is maturation, which is when the foliage starts to grow steadily.In this phenophase, the LAI varies with a low amplitude and can reach an annual peak.Subsequently, the foliage begins to enter into a recession, which is accompanied by an LAI decline.At a specific point in time, the foliage maturation corresponds to the end dates of the unfolding foliage.Obviously, this point is different from the vigorous period of unfolding foliage in the site observations.Thus, the large deviation between the retrieval results with the inversion algorithms and the observation data is shown in Figure 3.The phenophases derived from the inversion algorithms based on remote sensing data are also not consistent with the observed phenophases due to the different spatial scales of the data.The pixels of the satellite-derived phenophase primarily represent the compositive information of multiple vegetation types [2].However, ground observations on plant phenology are local observations covering only a few selected plant species [6].Thus, the difference between the phenological data is obvious.Furthermore, the area of the pixel based on the AVHRR LAI datasets is approximately 8 km × 8 km (a pixel of the MODIS LAI dataset is approximately 1 km × 1 km), which is considerably larger than the area of the site.This inconsistency may lead to the deviation between the retrieval results and the observation data, especially during senescence.Future in-depth research on reducing this inconsistency is urgently needed.

Shortage of Validation Data
Only three phenology sites were available in this study to validate the retrieval results for the whole phenophases from 1982 to 1988.Another ten observation sites were used to validate the retrieval results for the green-up phenophase, which included data ranging from 2000 to 2009.Clearly, the validation using the phenology observation data of minor sites was not very effective.In this study, we evaluated the inversion algorithms only at the single point scale, not at the regional scale, as the number of phenology observation sites is insufficient (less than 10).The shortage of observations increases difficulties greatly in the up-scaling process of phenophases from points to the region.Due to lack of reference maps of regional phenophases, validating the inversion algorithms at the regional scale is very difficult.Thus, enriching and perfecting the validation data in future research, including other geographical regions and biomes, is necessary to comprehensively evaluate the inversion algorithms.

Uncertainties of the Retrieval Results
Compared with the observation data of the phenology sites, the retrieval results using the inversion algorithms based on the remote sensing data revealed large uncertainties, which were related to the inversion algorithms and remote sensing datasets.The differences among the retrieval results using the inversion algorithms were analyzed.The differences in the retrieval results based on the various remote sensing datasets are shown in Figure 5. Three inversion algorithms show the retrieval results based on the improved MODIS LAI data were earlier than those based on the AVHRR LAI data.The bias statistics suggests the retrieved results based on the AVHRR LAI data are more reasonable than those based on the improved MODIS LAI datasets.So far, a comparison between the AVHRR LAI-3g and the improved MOIDS LAI has not been found in the literature, although many studies have compared and analyzed the NDVI datasets between AVHRR and MODIS, as well as their derivative growing season metrics [11,16,53].The strengths and weaknesses of the retrieval results based on the two LAI datasets were still unclear [11,54].Moreover, in the case of lacking the ground validation data, the comprehensive evaluation of the inversion algorithms is very difficult.Of the inversion algorithms, no single one can be accepted universally [22].Next, the temporally and spatially detailed historical datasets from both satellite and in situ observations are required to validate and correct the retrieval results in future plant phenology studies.
Additionally, the uncertainties of the retrieval results were related to the daily LAI derived from the data acquisition procedure.In this study, we used the cubic spline interpolation method to scale the LAI temporal frequency down from 15-day or 8-day to 1-day.The differences between the interpolated results and the real values always influence the retrieved results more or less even though this temporal interpolation method performs well.The uncertainties of the retrieval results were also related to the fitted LAI curve, although the S-curve function performed well in fitting the LAI curve.The methodology strongly relies on reasonable partitions of the LAI data.In this study, based on the previous published work in Shandong Province, the transition time separating the double seasonal crops was uniformly set to the 181st day of the year.This static threshold clearly does not apply to all stations, especially to the stations where agricultural cultivation patterns change frequently.Therefore, the unreasonable partition time of the LAI data in some stations leads to the unreliability of the fitted LAI curve and inaccuracy of the retrieval results.In view of the growing features of the seasonal vegetation, the crest number of the vegetation growth cycle could be confirmed by the fitted mathematical function [40].It may be more reasonable to divide the LAI data intervals with the time that corresponds to the trough between the crests of the vegetation growing curve.However, this partition time must still consult the farming work schedule in a given local area.Therefore, dividing the partition time of the LAI data dynamically must be further researched to enable the SC algorithm to be applied to a wider area.

Conclusions
Based on LAI datasets, we used the S-curve (SC) function integrating the asymmetric Gaussian (AG) function and the logistic (LC) function to fit the LAI curve.We then applied the resulting asymptotic lines and the curvature extrema to extract vegetation phenophases.The following conclusions were drawn by comparing with phenology observation data: (1) The SC algorithm has better performance for deriving the vegetation phenophases, especially for the germination and dormancy.The averaged absolute biases of germination are 2, 28 and 39 days for the SC, LC and AG algorithms, respectively; while the averaged absolute biases of dormancy are 5, 13 and 16 days for the SC, LC and AG algorithms, respectively.The retrieval results of the SC algorithm begin later than the other two algorithms during germination and green-up but start earlier during defoliation and dormancy.The discrepancy among the retrieval results is reflected in the fitted LAI curve.The SC function fits the LAI curve better than the other two functions.
The averaged RMSE values of fitting LAI are 0.116, 0.149 and 0.244 m 2 •m −2 for the SC, LC and AG functions, respectively.The averaged indices of agreement (IA) of fitting LAI are 0.983, 0.973, 0.892 for the SC, LC and AG functions, respectively.The phenology recognition rate (PRR) of the SC algorithm is obviously higher than the other two algorithms.The averaged PRRs are 1.00, 0.62, 0.91 for the SC, LC and AG algorithms, respectively.(2) Regardless of the LAI data used, the retrieved results using theSC algorithm are reliable and better than the other two algorithms.The averaged absolute biases of green-up using the SC algorithm are 9 and 12 days for the AVHRR LAI and improved MODIS LAI datasets, respectively.The green-up dates derived by the LC algorithm are earlier than the SC algorithm, but later than the AG algorithms.The SC and LC algorithms have the highest and lowest accuracy in the PRR (0.99 and 0.75), respectively.The SC algorithm has the advantage of deriving the vegetation phenology across time and space.(3) Besides the algorithms, the temporal-spatial resolution and quality of the LAI data used also determine the accuracy of the retrieved vegetation phenophases: low resolution/quality, higher uncertainty.For all three inversion algorithms, the retrieved results based on the AVHRR LAI data are later than those based on the improved MODIS LAI.The bias statistics analysis, i.e., the averaged absolute biases of green-up are 18 and 26 days for the AVHRR LAI and improved MODIS LAI data, respectively, shows that the retrieved results based on the AVHRR LAI data are more reasonable than those based on the improved MODIS LAI data.In addition, the mixing pixels of the LAI datasets and changes in the agricultural cultivation patterns influence the uncertainty as well.

Figure 2 .
Figure 2. Calculation schematic of vegetation phenophases (site position: 36.7°N,119.1°E; year: 2005).The LAI data are processed with the cubic spline interpolation.y l (y r ) is the fitted value of the left (right) S-curve function, and k l (k r ) is the corresponding curvature.

Figure 4 .
Figure 4. Comparison of fitted LAI using the S-curve function, logistic function and asymmetric Gaussian function for the year 1986 at the (a) Liaocheng, (b) Taian, and (c) Weifang site, respectively.The black dash line represents the transition time.

Figure 5 .
Figure 5.Comparison of green-up dates derived using the three inversion algorithms based on the AVHRR LAI data and improved MODIS LAI data at the (a) Heze, (b) Huimin, (c) Liaocheng, (d) Jining, (e) Taian, (f) Zibo, (g) Linyi, (h) Weifang, (i) Wendeng, and (j) Laiyang site, respectively.The superscript numbers 1 and 2 refer to the results based on the AVHRR LAI and the improved MODIS LAI, respectively.The superscript number 3 refers the dates were between early green-up and late green-up.

Table 1 .
Descriptions of phenology observation sites.

Table 2 .
Bias comparison of phenophase retrieval results (in days) using the three inversion algorithms.

Table 3 .
Comparison of phenological recognition rates (PRRs), root mean square error (RMSE) and index of agreement (IA) for fitting LAI (in m 2 •m −m ) among using the three inversion algorithms.