ANALYZING THE VELOCITY OF VEGETATION PHENOLOGY OVER THE TIBETAN PLATEAU USING GIMMS NDVI 3 g DATA

Global environmental change is rapidly altering the dynamics of terrestrial vegetation, and phenology is a classic proxy to detect the response of vegetation to the changes. On the Tibetan Plateau, the earlier spring and delayed autumn vegetation phenology is widely reported. Remotely sensed NDVI can serve as a good data source for vegetation phenology study. Here GIMMS NDVI3g data was used to detect vegetation phenology status on the Tibetan Plateau. The spatial and temporal gradients are combined to depict the velocity of vegetation expanding process. This velocity index represents the instantaneous local velocity along the Earth’s surface needed to maintain constant vegetation condition. This study found that NDVI velocity show a complex spatial pattern. A considerable number of regions display a later starting of growing season (SOS) and earlier end of growing season (EOS) reflected by the velocity change, particularly in the central part of the plateau. Nearly 74% vegetation experienced a shortened growing season length. Totally, the magnitude of the phenology velocity is at a small level that reveals there is not a significant variation of vegetation phenology under the climate change context. * Corresponding author


INTRODUCTION 1.1 Background
Under the global change context, how the spatial displacement of vegetation phenology changes correspond has attracted lots of attentions from the study filed (Elmore, 2016;Zhu, 2016).The current vegetation extent will be rearranged over the compact districts with climate change and human activities.The Tibetan Plateau region in China is a hotspot for both government and scholars are focusing on its environmental change and herdsman lives in the decades.As the third polar on the earth, Tibetan Plateau is sensitive to climate changes and human activities.Studies using satellite-derived data or experiment data, have pointed out it has experienced a continuous greening trend (Zhang, 2013;Shen, 2015).These research commonly applied spatial or temporal pattern analysis on the vegetation or NDVI-derived phenology.The concept of the velocity of change offers the opportunity to directly compare the ongoing change in the spatial patterns of vegetation phenology (Loarie, 2009).In this study, we mapped the velocity of change in vegetation phenology, which was calculated using the satellite-derived long-term NDVI.

GIMMS NVVI3g and Vegetation Phenology Extraction
The remotely sensed NDVI3g data used in this study was derived from NOAA's Advanced Very High Resolution Radiometer (AVHRR) sensor, which was produced in the framework of the Global Inventory Monitoring and Modeling System (GIMMS) project.This NDVI data has a spatial resolution of 8km under the WGS84 coordinate system and a temporal resolution of 16 days, covering the period of 1982-2012.More detailed information for data acquiring and processing is linked to references (Tucker, 2005).The regions where the long term mean NDVI is less than 0.05 are seemed as non-vegetation covered land and excluded in this study.On the Tibetan Plateau, the NDVI value ranges from 0.05 to 0.8 (Figure 1).The spatially averaged NDVI value for the multiyears mean over the study region is 0.36.
The 16-days repeated NDVI is constructed to a long term time series.Then, for each pixel, this time series is smoothed and interpolated with spline method.Based on the fitted growth curve, we calculated from time series annual metrics of vegetation phenology SOS (start of growing season), EOS (end of growing season), LOS (length of growing season).The White method that is a threshold-based method was applied to extract these phenology events (White, 1997).In this method, phenology metrics is calculated by scaling annual NDVI cycles between 0 and 1 (Equation 1).Vegetation growing season onset is observed at the day of year when NDVIratio exceeds 0.5.autumn offset is observed when NDVIratio falls below 0.5. (1) Where NDVImax and NDVImin is respectively the maximum and minimum NDVI in one year.After this scaled procedure, the observation NDVI values are scaled to the range 0-1, termed as NDVIratio.

Phenology Velocity
Both spatial and temporal pattern of vegetation phenology dynamic are the research key points depicted in the literature.It is a need to combine the spatial and temporal information together to characterize the vegetation phenology dynamic.Satellite based NDVI data is commonly used in the land dynamic process and climate changes studies from regional to global scale.The velocity concept was first developed for climate impact research to compare the displacement rate of a forcing variable (for example, temperature) with that of an affected variable (for example, the occurrence of a given species) by converting them into the same units, such as km/yr (Huang, 2017).In this case, the velocity of phenology change indicates a displacement in the isolines of this variable due to climate change for multiple pixels in a certain region.Here, phenology velocity index of change (km/yr) was derived from spatial gradient (doy/km) and the temporal linear trend of phenology for 31 years (doy/yr).This index is inspired by the global climate change velocity for vegetation activity [6]. (2) (3) Where the temporal trend is derived from the slope of linear regression between phenology metrics time series and the corresponding years for each pixel.The pixel with trend that pass the 0.05 significant level test is used for the velocity calculation, otherwise is excluded.
The spatial gradient for each pixel in NDVI and phenology image is calculated as Equation.4.We chose a template for 3*3 to calculate the gradient for the central pixel, which will be impacted by its eight neighbours. (4) Where is the spatial gradient for pixel with the x, y axis coordinate (i, j).dx and dy are horizontal direction and vertical direction convolution kernel, respectively, determined as following Equation.5.
(5) dx kernel is used to detect horizontal gradient, and dy kernel is used to detect vertical gradient.The central pixel gradient is affected by its neighbors, but varying among their distances.The closer neighbour is assigned to a bigger weight.For example, in dx kernel matrix, central pixel's horizontal neighbours have a weight of 2 or -2, and neighbours on diagonal line is weighted by 1 or -1.It is the same as the dy matrix.Following procedure computed the gradient for central pixel through NDVI and phenology values multiplying horizontal and vertical kernels (Equation 4).

Spatial Pattern of NDVI and Phenology Trends
According to mulitple years mean of NDVI, nearly half of the vegetated area on the Tibetan Plateau has a relatively high vegetation cover, see Figure 1.The spatially averaged NDVI is 0.26.The northern arid area is mostly covered by bare land.The vegetation is mainly distributed in the eastern and southern regions where have a lower elevation and better natural environment.With regard to the long term trend of NDVI and phenology evolution, their pixel level temporal gradient was computed as the slope of a linear model fit through each year of the time period of interest, at the 5% significance level (Figure 2).The spatial pattern of NDVI is similar to that of mulitple years mean NDVI.The majority of the vegetated area show an increasing trend of NDVI.But there are some regions experienced a vegetation degrading process, particularly in the southern Tibetan Plateau (Figure 2a).For the phenology metrics trend analyses, the pixel value that is failing to extract phenological parameters is set as nodata in Figure 2c-2d.For SOS, there are a considerable number of pixels showing a inceasing trend in the southern and central part.This phenomenon indicates that plants greenup dates are delayed in the springtime.Whereas, plants andvanced greenup dates are reflected by the deceasing slope of the SOS trends in the eastern and western part of the plateau (Figure 2b).For EOS, positive trend means the delayed cessation in the fall.While negative EOS slope reflects the advanced offset of growing season.There are more than 75% of EOS trend pixels showing a negative value.Based on the SOS and EOS dates, we calculated the LOS length, in which positive slope indicates a prolonged growsing season period.In the northwestern region and part of the southern region, LOS display a shortened condition (red color in Figure 2d, about 67% in the valid pixels).

Spatial Pattern of Phenology Velocity
Based on the spatial and temporal gradients, we calculated the velocity of vegetation phenology for the Tibetan Plateau (Figure 3).The spatially averaged NDVI velocity for the study area is 0.014km/yr, with the standard deviation of 0.16km/yr.Compared with the spatial pattern of multi-years mean of NDVI, the spatial distribution velocity covers a less area because the phenology extraction method filters out the areas where vegetation growing curve is not fitted well.For NDVI velocity, the majority of the area show a positive value indicating the greening trend (36% in valid pixels), but the velocity magnitude is small.NDVI velocity show a highly spatial heterogeneity pattern.
The average velocity of change in SOS, EOS and LOS over the study area is 0.023± 0.17km/yr, -0.044 km/yr, -0.078 km/yr, respectively.SOS display a positive velocity of change in the southwestern part of Tibetan Plateau (36%), and some negative velocities in the eastern part.Compared with SOS, a large part of study areas (74.3% in valid pixels) show a negative EOS velocity, revealing the earlier vegetation offset in the autumn.As a result of the combination of SOS and EOS, we obtained the spatial distribution of LOS velocity (Figure 3d).LOS spatial pattern is similar to that of EOS, but with more area showing a negative velocity of change (nearly 74%).

CONCLUSION
In summary, in this study we have applied the concept of velocity to remotely sensed vegetation phenology fields which combines spatial and temporal gradients information to detect the vegetation phenology dynamic on the Tibetan Plateau.Although there is an obvious expansion of the phenology dynamic over the last three decades, the activity magnitude is not very large, representing by the low velocities of SOS, EOS and EOS change.There is a majority number of pixels displaying a negative velocity of phenology change.About 72.7% of the valid EOS and 74.3% of the LOS show a negative velocity of phenology change.SOS and NDVI has a smaller value of 36% and 36.13%respectively.Although totally the Tibetan Plateau displays a greening trend in the literatures, it not so significant in the detection of velocity changes.

Figure 1 .
Figure 1.Multi-years mean of NDVI over Tibetan Plateau

Figure 3 .
Figure 3. Velocity for northeastern China.The velocity unit is km/yr.