Assessing Performance of NDVI and NDVI3g in Monitoring Leaf Unfolding Dates of the Deciduous Broadleaf Forest in

Using estimated leaf unfolding data and two types of Normalized Difference Vegetation Index (NDVI and NDVI3g) data generated from the Advanced Very High Resolution Radiometer (AVHRR) in the deciduous broadleaf forest of northern China during 1986 to 2006, we analyzed spatial, temporal and spatiotemporal relationships and differences between ground-based growing season beginning (BGS) and NDVI (NDVI3g)-retrieved start of season (SOS and SOS3g), and compared effectiveness of NDVI and NDVI3g in monitoring BGS. Results show that the spatial series of SOS (SOS3g) correlates positively with the spatial series of BGS at all pixels in each year (P < 0.001). Meanwhile, the time series of SOS (SOS3g) correlates positively with the time series of BGS at more than 65% of all pixels during the study period (P < 0.05). Furthermore, when pooling SOS (SOS3g) time series and BGS time series from all pixels, a significant positive correlation (P < 0.001) was also detectable between the spatiotemporal series of SOS (SOS3g) and BGS. In addition, the spatial, temporal and spatiotemporal differences between SOS (SOS3g) and BGS are at acceptable levels overall. Generally speaking, SOS3g is more consistent and accurate than SOS in capturing BGS, which suggests that NDVI3g data might be more sensitive than NDVI data in monitoring vegetation leaf unfolding.


Introduction
Monitoring the vegetation growing season is crucial not only for assessing ecosystem responses to climate change [1,2] but also for identifying the carbon-uptake period [3][4][5][6][7][8][9] and examining the seasonal exchanges of water and energy between land surface and atmosphere [10,11].So far, there are two types of approaches for detecting the vegetation growing season.The conventional phenology approach has focused on observing and simulating the ground-based growing season at spatially discrete sites using the phenological data of individual plants [12][13][14][15][16].This kind of approach identifies the individual plant growing season at local scales but may not represent the growing season of plant communities at regional scales.By contrast, the land surface phenology approach determines the satellite-derived growing season retrieved from various vegetation index data, such as the Normalized Difference Vegetation Index data from the Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) [17][18][19][20][21][22][23][24][25][26][27][28] over broad and continuous geographic coverage.Because metrics and thresholds of vegetation indices may not directly correspond to conventional, ground-based phenological events, but rather provide indicators of vegetation dynamics [19], a detailed comparison of these satellite measures with ground-based phenological events is needed [29][30][31][32][33].However, validation of remote sensing retrievals is an important but difficult task because of spatial scale differences, temporal resolution of the observations, and quality of the field measurements [32,33].In recent years, some studies have been carried out to compare satellite-derived onset and offset of greenness with surface phenological stages of individual plants at selected locations, and found relatively consistent spatiotemporal variations between ground-based and satellite-derived growing season beginning and end dates [32][33][34][35].Nevertheless, this kind of validation of land surface phenology might be limited by a lack of in situ datasets collected across ecosystems at spatial scales commensurate with satellite data [31,33,[36][37][38][39][40][41][42].Substantial validation is therefore possible by upscaling intensive field phenology measurements to vegetation community and landscape scales [42].On the basis of above considerations, we established gridded spatial datasets of growing season beginning dates year-by-year by extrapolating leaf unfolding (LU) dates of indicator tree species from the stations of model fitting to all pixels within the research region using station-specific spatial phenology models [16] and gridded daily temperature data in each year.In this way, we obtained a ground-based phenology dataset at spatial scales commensurate with the satellite-derived phenology dataset, so that the pixel to pixel comparison can be implemented.Both AVHRR Normalized Difference Vegetation Index (NDVI) data and AVHRR third generation Normalized Difference Vegetation Index (NDVI3g) data were used for this research.The main objective was to compare performance of the satellite-derived start of season (SOS) from AVHRR NDVI and NDVI3g in monitoring ground-based growing season beginning (BGS) based on spatial and temporal correlation analysis and error estimates.

Study Area and Plant Species Selection
Many studies have shown that performance of land surface phenology in monitoring ground vegetation phenology depends significantly on selected vegetation types [26,32,[43][44][45][46].Since seasonal feature change of the temperate deciduous broadleaf forest (TDBF) is remarkable in the Northern Hemisphere, we selected the TDBF areas in northern China to analyze the relationship between the ground observed phenological stage and remotely sensed deciduous canopy development.In order to obtain real and precise vegetation information, we used land cover dataset to determine spatial ranges of the TDBF as the study areas (Figure 1).The TDBF is a dominant forest vegetation type in northern China, which is distributed mainly in mountainous areas (Figure 1).The plant community is composed of heliophilous broadleaf trees and shrubs with remarkable seasonality, including Quercus, Betula, Carpinus, Ulmus, Celtis, Acer, Populus, Malus, etc. [47].In order to establish an in situ dataset of spring phenology across the TDBF, we selected four dominant tree species to represent the local plant community.The selected tree species include Salix matsudana (Hankow Willow), Populus simonii (Simon Poplar), Ulmus pumila (Siberian Elm) and Prunus armeniaca (Common Apricot) that are all native deciduous trees and grow widely in temperate northern China.In this study, we defined mean LU beginning dates of the four dominant tree species as ground-based BGS of local plant community.

Phenological Data
The LU data during 1986 to 2006 were acquired from the China Meteorological Administration (CMA).The CMA network is the largest phenological observation system in China and came into operation in 1980 [48].Observation contents include phenophases from 28 common species of woody plants, one common species of herbaceous plant, and 11 common species of animals [49].Among the observed woody plants, the number of stations and the recorded data for the four tree species are much more than other species.According the observation criteria [49], LU beginning was identified when a few leaves are fully spreading in spring.It should be noted that because time series lengths of LU beginning data are not identical among stations and species, the number of stations for LU modeling is different for the four tree species, namely, 77 stations for Salix matsudana, 61 stations for Populus simonii, 72 stations for Ulmus pumila, and 40 stations for Prunus armeniaca, which are located in temperate northern China (Figure 2).

Climate Data
Daily mean air temperature data at 343 stations in temperate northern China from 1986 to 2006 were acquired from the China Meteorological Data Sharing Service System (http://cdc.cma.gov.cn/).In order to establish in situ dataset of LU beginning over a continuous geographic coverage, we used climate data interpolation package ANUSPLIN 4.2 [50] and Digital Elevation Model (DEM) data derived from the US Geological Survey to interpolate the daily mean air temperature into 8 km × 8 km grids over the deciduous broadleaf forest areas.Based on the gridded daily mean temperature dataset and daily temperature-based spatial phenology models [16], we estimated LU beginning dates of the four tree species at each grid and in each year.

Satellite Data and Data Processing
Two types of NDVI datasets generated from AVHRR sensor by Global Inventory Monitoring and Modeling System (GMMIS) [51] were used in this study.The available NDVI data are from July 1981 to December 2006 with the spatial resolution of 8 km and temporal interval of 15 days.So far, the AVHRR NDVI is the most widely used data source for monitoring vegetation dynamics over large geographical coverage during the past two decades [22,26,28,52,53].The AVHRR NDVI3g dataset was also developed under the framework of GIMMS project.The spatial resolution and temporal interval of the AVHRR NDVI3g dataset are similar to AVHRR NDVI but the time series length is prolonged to December 2011.AVHRR NDVI3g was carefully assembled from different AVHRR sensors and accounted for several deleterious effects, such as calibration loss, orbit drift, volcanic eruption.
Land cover data used in this study is the Global Land Cover Classification collection (GLCC) generated by the University of Maryland [54].The GLCC product utilized the taxa method proposed by the International Geosphere-Biosphere Program (IGBP) and was validated by Landsat Thematic Mapper (TM) images [55].Imagery acquired from the AVHRR satellites between 1981 and 1994 was analyzed to distinguish fourteen land cover classes.Since this product is available at 8 km pixel resolutions, it is suitable for defining land cover types in our study.
With regard to data processing, first, we used the Savitzky-Golay filter to smooth the original NDVI time series so that abnormal high and low values could be eliminated [56].Then, we selected midpoint model [22] to extract the satellite-derived start of season at each pixel and in each year.Here, the start of season derived from AVHRR NDVI dataset was denoted as SOS, while the start of season derived from AVHRR NDVI3g dataset was denoted as SOS3g.According to ranges of observed LU dates of dominant tree species in northern China, we eliminated the abnormal SOS and SOS3g data outside the period of 50-180 day of year (DOY).

Daily Temperature-Based Spatial Phenology Model
The basic hypothesis is that the station-to-station variation of a phenological event occurrence date over an area is mainly influenced by station-to-station variation of daily mean temperature within a particular length period (LP) of days during and before its occurrence over the area.In order to determine the LP during which the station-to-station variation of daily mean temperature affects station-to-station variation of the phenological occurrence date most remarkably in a year, we first calculated the number of days between the earliest and latest date in spatial series of the phenological occurrence dates across the stations in the year, calling this the basic LP (bLP).Then, we computed the daily mean temperature spatial series at the stations in the year during the bLP plus a moving LP (mLP) prior to the earliest date in spatial series of the phenological occurrence date by step length of one day, namely, during bLP + 1 day, bLP + 2 days, bLP + 3 days, etc.Thus, the LP is defined as follows: Further, we calculated correlation coefficients between spatial series of the phenological occurrence date and spatial series of daily mean temperature during different LPs (bLP + 1 day, bLP + 2 days, bLP + 3 days, etc.) at the stations in the year.Finally, we obtained the optimum LP with the largest correlation coefficient between spatial series of the phenological occurrence date and spatial series of daily mean temperature in the year [16].
Based on the species-specific daily mean temperature within the optimum LP in each year, we established daily mean temperature-based spatial LU regression models for the four tree species and in each year.External validation showed that the daily mean temperature-based spatial LU regression models have a strong spatial extrapolation capability to LU dates of the four tree species at external stations.Since these models allowed simulation and prediction of yearly spatial patterns of LU dates on the basis of yearly spatial patterns of LP temperatures, we substituted the species-specific LP temperatures at all 8 km × 8 km grids within the deciduous broadleaf forest areas in each year into the corresponding yearly spatial LU models of the four tree species, respectively, and obtained gridded ground-based LU dates of the four tree species in each year.Then we calculated mean LU dates of the four tree species at each grid and in each year as the ground-based LU dataset.Figure 3 shows procedures of daily mean temperature-based spatial LU simulation and extrapolation.

Error Estimate
In order to estimate the difference between satellite-derived SOS (SOS3g) and ground-based BGS, we used two criteria: the mean error (ME) between SOS (SOS3g) and BGS, and the mean absolute error (MAE) between SOS (SOS3g) and BGS.They were calculated by the following formula: where SOSi is the satellite-derived start of season (DOY) in year i or at pixel i; BGSi is the ground-based beginning of the growing season (DOY) in year i or at pixel i; n is the number of years for estimating the difference between SOS (SOS3g) and BGS time series or the number of pixels for estimating the difference between SOS (SOS3g) and BGS spatial series.

Spatial Relationship between Satellite-Derived SOS (SOS3g) and Ground-Based BGS
Generally speaking, mean BGS, SOS and SOS3g show a similar spatial pattern with an obviously latitudinal progression from south (earlier) to north (later) over 1986 to 2006.In addition, an approximately longitudinal progression can also be detected from east (earlier) to west (later) in the eastern part of northern China.The spatial differences in BGS, SOS and SOS3g were from 90 DOY (the end of March) to 170 DOY (early June) (Figure 4).In order to precisely assess spatial consistency between SOS (SOS3g) and BGS, we carried out a correlation analysis between spatial series of SOS (SOS3g) and BGS year-by-year.Results show that both SOS and SOS3g correlates positively with BGS in each year (P < 0.001).We noted that spatial correlation coefficients between SOS3g and BGS are larger than those between SOS and BGS in each year, which indicates that spatial patterns of AVHRR NDVI3g-derived SOS may be more consistent with spatial patterns of ground-based BGS (Table 1).
With regard to mean differences between satellite-derived SOS (SOS3g) and ground-based BGS at all pixels in each year, the MEs are between −2.8 days (1998) and 6.1 days (1990) for SOS and between −6.2 (2004) and 1.9 days (1990) for SOS3g.The absolute values of MEs are less than 5 days in 19 of 21 years (90%) for both SOS and SOS3g.Meanwhile, the MAEs are between 6.6 (1988) days and 11.0 days (2004) for SOS and between 6.0 days (1988) and 10.0 days (2004) for SOS3g.Overall, the MAEs for SOS3g are obviously smaller than the MAEs for SOS in all year except 2005.That is, the AVHRR NDVI3g-derived SOS may be more similar to the ground-based BGS than the AVHRR NDVI-derived SOS (Table 1).

Temporal and Spatiotemporal Relationship between Satellite-Derived SOS (SOS3g) and Ground-Based BGS
A significantly positive correlation (P < 0.05) was found between SOS and BGS at 65.8% of all pixels and between SOS3g and BGS at 67.3% of all pixels during 1986 to 2006.Nonsignificant correlation coefficients are mainly located in the south part of the Da Hinggan Ling of northeastern China (Figure 5(a,b)).Therefore, the temporal consistency between satellite-derived start of season and ground-based beginning of the growing season existed in most parts of northern China and was slightly larger between SOS3g and BGS than between SOS and BGS.Considering mean differences between SOS (SOS3g) and BGS at each pixel during 1986 to 2006, the mean MEs are between +10 days and −10 days at 74.3% of all pixels for SOS and 79.4% of all pixels for SOS3g.Spatial patterns show that positive MEs were detected mostly in south parts of northern China, whereas negative MEs were evident mostly in north parts.Namely, SOS and SOS3g were generally later than BGS in south parts but earlier in north parts (Figure 5(c,d)).Moreover, the mean MAEs are less than 10 days at 70.6% of all pixels between SOS and BGS and 76% of all pixels between SOS3g and BGS.The lowest values of mean MAEs appeared in the Changbai Shan (Figure 5(e,f)).In general, the mean differences between SOS3g and BGS are slightly smaller than those between SOS and BGS.Furthermore, the spatiotemporal analysis shows that both SOS and SOS3g correlate positively with BGS (P < 0.001), namely, the earlier the BGS at a pixel in a year, the earlier the SOS and SOS3g at the pixel in the year.Besides, all observation pairs are concentrated near by the 1:1 diagonal line, which indicates that the spatiotemporal variations of SOS (SOS3g) and BGS have consistent characteristics.Again, the spatiotemporal correlation coefficient between SOS3g and BGS is larger than that between SOS and BGS, while the MAE between SOS3g and BGS are smaller than that between SOS and BGS (Figure 6).

Discussion
In this study, gridded LU estimation data of the four tree species were used to validate AVHRR NDVI (AVHRR NDVI3g)-retrieved start of season in the deciduous broadleaf forest of northern China during 1986 to 2006.Comparing our study with similar validation studies over the past 20 years worldwide, spatiotemporal correlation coefficients and MEs between SOS (SOS3g) and BGS in northern China fall in the range of those in other regions (Table 2).Nevertheless, our study validated the satellite-derived start of season over a broader geographical coverage and longer time series than other studies so far.
Comparing spatial, temporal and spatiotemporal patterns of BGS with those of SOS (SOS3g), we found that the SOS3g is more similar to the BGS than the SOS (Table 1, Figures 5 and 6).In order to assess sensitivity of NDVI and NDVI3g products in monitoring vegetation LU, we calculated regional mean AVHRR NDVI and AVHRR NDVI3g values for each 15-day period during the 1986-2006 period.Figure 7(a) shows that the interannual fluctuations of both datasets are consistent but the interannual amplitudes of them are obviously different, especially, the peak values of NDVI3g are significantly higher than those of NDVI (Figure 7(b)).This implies that NDVI3g data might be more sensitive than NDVI data in monitoring vegetation LU.Therefore, the spatial, temporal and spatiotemporal correlation coefficients between SOS3g and BGS are larger than those between SOS and BGS.As mentioned above, spatial series of SOS (SOS3g) correlate positively with spatial series of BGS in each year (P < 0.001).However, spatial correlation coefficients between SOS (SOS3g) and BGS indicate an obviously interannual variation (Table 1), which may be associated with interannual variation in the representativeness of mean LU date of the four tree species to green-up date of the local plant communities.We assume that the range (days) between the earliest and latest LU dates of the four tree species in a year can reflect the time period of green-up process of the local plant communities in the year to a certain extent.Therefore, a larger range implies a prolonged green-up process of the local plant communities, while a smaller range means a shortened green-up process of the local plant communities.Therefore, the uncertainty of mean LU date of the four tree species in representing green-up date of the local plant communities might be enhanced in the years with a larger range and reduced in the years with a smaller range.According to above consideration, we may expect that the smaller the mean range between the earliest and latest LU dates of the four tree species at all pixels in a year, the better the representativeness of mean LU date of the four tree species to green-up date of the local plant communities in the year, which may induce a larger spatial correlation coefficient between SOS (SOS3g) and BGS in the year.Statistical analysis confirmed the above assumption and expectation, namely, the yearly spatial correlation coefficients between SOS (SOS3g) and BGS correlate negatively with the yearly mean range between the earliest and latest LU dates of the four tree species at all pixels from 1986 to 2006.Therefore, the representativeness of mean LU date of dominant tree species to green-up date of the local plant communities is an important influencing factor on effectiveness of validating satellite-derived start of season (Figure 8).

Conclusions
This study validated AVHRR NDVI (AVHRR NDVI3g)-retrieved start of season by means of gridded mean LU dates of four tree species in the deciduous broadleaf forest of northern China during 1986 to 2006.The main conclusions are as follows: 1.The spatial series of SOS (SOS3g) correlates positively with the spatial series of BGS across the research region in each year (P < 0.001).Meanwhile, the time series of SOS (SOS3g) correlates positively with the time series of BGS at more than 65% of all pixels during the study period (P < 0.05).When pooling SOS (SOS3g) time series and BGS time series from all pixels, a significant positive correlation (P < 0.001) was also detected between the spatiotemporal series of SOS (SOS3g) and BGS.In addition, the spatial, temporal and spatiotemporal differences between SOS (SOS3g) and BGS are at acceptable levels overall.Therefore, the NDVI (NDVI3g)-retrieved start of season can effectively capture spatial patterns and temporal variations of the ground-based growing season beginning.2. The NDVI3g-retrieved start of season is generally more consistent and accurate than the NDVI-retrieved start of season in capturing the ground-based growing season beginning.This finding suggests that NDVI3g data might be more sensitive than NDVI data in monitoring vegetation LU. 3. The representativeness of mean LU date of dominant tree species to green-up date of the local plant communities is an important influencing factor on effectiveness of validating satellite-derived start of season.

Figure 1 .
Figure 1.Distribution of the temperate deciduous broadleaf forest areas in northern China.

Figure 2 .
Figure 2. Distribution of phenological stations and time series length of leaf unfolding data for the four tree species.

Figure 3 .
Figure 3. Flow-diagram of daily mean temperature-based spatial LU simulation and extrapolation.

Figure 4 .
Figure 4. Spatial pattern of mean BGS (growing season beginning), SOS (start of season) and SOS3g over 1986 to 2006 in the deciduous broadleaf forest of northern China.(a) BGS; (b) SOS; (c) SOS3g.

Figure 5 .
Figure 5. Spatial patterns of temporal correlation coefficients and mean differences between SOS (SOS3g) and BGS at each pixel from 1986 to 2006.(a) correlation coefficients between SOS and BGS; (b) correlation coefficients between SOS3g and BGS; (c) ME between SOS and BGS; (d) ME between SOS3g and BGS; (e) MAE between SOS and BGS; (f) MAE between SOS3g and BGS.

Figure 6 .
Figure 6.Spatiotemporal correlation coefficients and differences between SOS (SOS3g) and BGS at each pixel in each year.(a) between SOS and BGS; (b) between SOS3g and BGS.

Figure 7 .
Figure 7. (a) Interannual variations of regional mean AVHRR NDVI and AVHRR NDVI3g and (b) relationships between regional mean AVHRR NDVI and AVHRR NDVI3g from 1986 to 2006.

Figure 8 .
Figure 8. Relationship between SOS (SOS3g)-BGS spatial correlation coefficient and the mean range of leaf unfolding dates of the four tree species at all pixels from 1986 to 2006.(a) between SOS-BGS spatial correlation coefficients and the mean range of leaf unfolding dates; (b) between SOS3g-BGS spatial correlation coefficients and the mean range of leaf unfolding dates.

Table 1 .
Spatial correlation coefficients (r), mean error (ME) and mean absolute error (MAE) between SOS (SOS3g) and BGS at all pixels in each year (n = 3,994).

Table 2 .
Comparison in spatiotemporal validation of satellite-derived spring phenology based on different Vegetation Indices and conventional phenological data.