An Exploration of Terrain Effects on Land Surface Phenology across the Qinghai–Tibet Plateau Using Landsat ETM+ and OLI Data

: Detecting spatial patterns of land surface phenology (LSP) with high spatial and temporal resolutions is crucial for accurately estimating phenological response and feedback to climate change and biogeochemical cycles. Numerous studies have revealed LSP across the Qinghai–Tibet Plateau (QTP) using a variety of coarse-resolution satellite data. However, detailed phenological spatial patterns along with changes of mountain topography remain poorly understood, which greatly limits efforts to predict the impacts of climate change on vegetation growth and ecosystem productivity in complex terrain regions. Combining Landsat 7 ETM+ (Enhanced Thematic Mapper Plus) and Landsat 8 OLI (Operational Land Imager) observations in overlapping zones of adjacent images, this study detected Normalized Difference Vegetation Index (NDVI)-based LSP metrics at a spatial resolution of 30 m, and explored how LSP varied with topographic factors along a 2600 km belt transect of the central QTP. The results show that the greenup onset date showed a delayed tendency with the increase of elevation at a mean rate of 1.52 days/100 m, while the dormancy onset date indicated an advanced tendency at a mean rate of − 0.59 days/100 m. In general, greenup onset date was later but dormancy onset date was earlier on shaded slopes than on sunlit slopes in the meadow area. By contrast, greenup onset date did not signiﬁcantly depend on aspect in the steppe area, while dormancy onset date indicated a similar response to aspect in the steppe area with that in the meadow area. With regard to the effect of slope on vegetation phenology in the meadow area, greenup onset date was signiﬁcantly delayed but dormancy onset date signiﬁcantly advanced with the increase of slope on both north and south slopes. In the steppe area, however, the inﬂuence pattern of slope on vegetation phenology was the opposite. Essentially, effects of topographical parameters on LSP were controlled by temperature and moisture combinations in complex terrain.


Introduction
Phenology is the study of recurring plant and animal life cycle stages, especially their timing and relationships with weather and climate [1]. Plant phenology focuses on the timing of annually recurring plant growth and reproductive phenomena, as well as the drivers of these events associated with endogenous and exogenous forces [2]. Plant phenological variation is highly sensitive to climate change and has strong feedbacks to the climate system [3]. Since the Qinghai-Tibet Plateau (QTP) has experienced a rapid warming that is twice as high as the global average level in recent decades [4], climate change and its impacts on alpine vegetation phenology have been the source of much concern [5][6][7][8][9]. To date, most land surface phenological studies on the QTP were based on Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Leaf Area Index (LAI) data from the Advanced Very High Resolution Radiometer (AVHRR) [5,6,[10][11][12][13][14][15][16][17], Moderate-Resolution Imaging Spectroradiometer (MODIS) [6][7][8]14,18], and Satellite Pour l'Observation de la Terre (SPOT) [6][7][8]14,[19][20][21]. The spatial resolutions of these data range normally from 500 m to 8 km. Thus, the extracted phenological metrics generally represent the composite of various vegetation types in a relatively large terrain unit and have been used to detect general phenological response to climate change at regional scales. By contrast, detailed spatial phenological patterns and their response to localized climatic factors remain poorly understood, which might lead to an inaccurate estimate of carbon budgets and ecosystem productivity on the QTP. To investigate spatial phenological responses to localized climatic factors, topographic features (such as slope, aspect, and elevation) at higher spatial resolution (such as 30 m) that control vegetation properties and growth conditions at local areas [22][23][24] can be used as a proxy for localized climatic factors. This is based on the fact that the spatial resolutions of gridded climate datasets currently available are normally coarser than 0.05 degrees [25], which is unable to represent the microclimate differentiation on the highly heterogeneous QTP.
Landsat sensors observe the land surface at a spatial resolution of 30 m with a repeat cycle of 16 days. However, Landsat observations with low cloud coverage are very limited within a year for most areas on the Earth [26]. Therefore, multiyear Landsat NDVI/EVI data have been combined to improve the temporal resolution for phenological detection [27][28][29]. This kind of phenological detection can capture multiyear mean vegetation dynamics at higher spatial resolution than that with satellite data from other sensors with more frequent repeat cycles. In addition, interannual variation of vegetation phenology could be derived by computing the annual anomaly in the timing (day of year) of the spring onset and autumn onset relative to the long-term average [29,30]. This approach can produce reasonable interannual phenological variation in North American temperate and boreal deciduous forests [28,29,31], and mixed and conifer-dominated forests [31]. Alternatively, the temporal resolution of Landsat data can also be improved by fusing Landsat streams with daily MODIS data for detailed analysis of vegetation phenology [32][33][34]. The basic assumption in current data fusion methods is that the vegetation growth rate is the same within a land cover type [35]. It should be noted that although these approaches may roughly retrieve multiyear mean and yearly land surface phenology from Landsat data, accurate vegetation phenology timing metrics should be derived from Landsat observations in the same year. Therefore, our research developed a new approach to estimate the grassland LSP using a single-year NDVI curve on the QTP.
The overall objectives of this study are to (1) extract vegetation phenology at each pixel along a transect across the middle of the QTP based on the Landsat observations in the same year; (2) quantify the impacts of topographic features on spatial heterogeneity of vegetation phenological metrics; and (3) explain the climatic cues of terrain effects. To achieve these goals, we combined Landsat 7 ETM+ (Enhanced Thematic Mapper Plus) and Landsat 8 OLI (Operational Land Imager) at the overlapping zones of adjacent images across the middle of the QTP from 2013 to 2015. This combination could provide as many as 92 observations within a year, which is potentially able to produce fairly accurate phenological timings [36].

Study Area and Vegetation Types
The Qinghai-Tibet Plateau is located in western and southwestern China, covering approximately 2.3 million km 2  increases rapidly from about 2000 m in the east to more than 8000 m in the west with an average elevation exceeding 4000 m. Influenced by both elevation and Indian monsoon, the climate shows a thermal/moisture gradient from southeast to northwest. According to the climate regionalization scheme [37], the climate of the QTP can be divided into plateau's temperate humid region, plateau's temperate and subalpine subhumid region, plateau's temperate and subalpine semiarid region, and plateau's temperate and alpine arid region. The annual mean air temperature decreases from 15.5 • C in the southeast to −5.0 • C in the northwest and the annual mean total precipitation decreases from 1764 mm in the southeast to 16 mm in the northwest for the period of 1981 to 2011 [9]. The dominant vegetation types along the thermal/moisture gradient include forest, shrub, meadow, steppe, and desert [38].
In this study, we focused on steppe and meadow ( Figure 1). These two major vegetation types are distributed across the QTP except for southeastern forest/shrub areas and northwestern desert areas, and cover about 60% of the QTP. Specifically, meadow is mainly concentrated in the central, eastern, and southern areas with temperate and subalpine subhumid climate, while steppe appears mostly in the western and northern areas with temperate and subalpine semiarid and arid climate.
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 17 on the QTP increases rapidly from about 2000 m in the east to more than 8000 m in the west with an average elevation exceeding 4000 m. Influenced by both elevation and Indian monsoon, the climate shows a thermal/moisture gradient from southeast to northwest. According to the climate regionalization scheme [37], the climate of the QTP can be divided into plateau's temperate humid region, plateau's temperate and subalpine subhumid region, plateau's temperate and subalpine semiarid region, and plateau's temperate and alpine arid region. The annual mean air temperature decreases from 15.5 °C in the southeast to −5.0 °C in the northwest and the annual mean total precipitation decreases from 1764 mm in the southeast to 16 mm in the northwest for the period of 1981 to 2011 [9]. The dominant vegetation types along the thermal/moisture gradient include forest, shrub, meadow, steppe, and desert [38].
In this study, we focused on steppe and meadow ( Figure 1). These two major vegetation types are distributed across the QTP except for southeastern forest/shrub areas and northwestern desert areas, and cover about 60% of the QTP. Specifically, meadow is mainly concentrated in the central, eastern, and southern areas with temperate and subalpine subhumid climate, while steppe appears mostly in the western and northern areas with temperate and subalpine semiarid and arid climate.

Satellite Data and Data Processing
A total of 4135 Landsat ETM+ (Enhanced Thematic Mapper Plus) and OLI (Operational Land Imager) images from 2013 to 2015 were obtained from the United States Geological Survey (USGS) archive (http://earthexplorer.usgs.gov/). Based on these data, 17 Landsat overlapping zones of adjacent images were selected along a transect across the middle of the QTP following the elevational gradient from east to west ( Figure 1). Fewer OLI images are available in 2013 than in 2014 and 2015 because Landsat 8 launched in February 2013 ( Figure A1 in Appendix A). To reduce atmospheric effects, raw digital number values from each image were converted to surface reflectance using the atmosphere correction tool in the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [39,40]. Surface reflectance data from both Landsat 7 and 8 images were used to calculate the Normalized

Satellite Data and Data Processing
A total of 4135 Landsat ETM+ (Enhanced Thematic Mapper Plus) and OLI (Operational Land Imager) images from 2013 to 2015 were obtained from the United States Geological Survey (USGS) archive (http://earthexplorer.usgs.gov/). Based on these data, 17 Landsat overlapping zones of adjacent images were selected along a transect across the middle of the QTP following the elevational gradient from east to west ( raw digital number values from each image were converted to surface reflectance using the atmosphere correction tool in the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [39,40]. Surface reflectance data from both Landsat 7 and 8 images were used to calculate the Normalized Difference Vegetation Index (NDVI). In addition, we used the Fmask product [41,42] to remove cloud, cloud shadow, snow cover, and open water.
NDVI from Landsat 7 ETM+ was calibrated to be comparable to that from Landsat 8 OLI. A uniform transformation function between OLI NDVI and ETM+ NDVI values was developed based on datasets in the United States, which was assumed to be applicable for various regions and time periods [43]. However, we found that the transformation function varied considerably on the QTP. Thus, we established a set of temporally local functions using linear regressions between spatial OLI and ETM+ NDVI data pairs in temporally adjacent time phases in each zone. The temporally local functions were applied to calibrate Landsat 7 ETM+ NDVI.
The MODIS land surface temperature (LST) product (MOD11A2) with a temporal resolution of 8 days and spatial resolution of 1000 m from 2013 to 2015 was also acquired from the USGS archive (http://earthexplorer.usgs.gov/), and was then used to determine the winter period based on a threshold of 278 K. The maximum NDVI during the winter was taken as the background value for vegetation dormancy because temporal NDVI variation below 278 K reflects snow contamination [44]. To match the spatial resolution of Landsat observations, LST data were resampled into 30 m resolution using a nearest neighbor method.

Land Surface Phenology Extraction Using Landsat Data
The Hybrid Piecewise Logistic Model (HPLM) algorithm [44], which was developed by combining a piecewise logistic model [45] with greenness under stress [29,46], has been used to detect phenological metrics from AVHRR, MODIS, and VIIRS (Visible Infrared Imaging Radiometer Suite) data [47,48]. A sensitive analysis indicated that vegetation phenology can be estimated precisely from satellite data with temporal resolutions of 6-16 days [36]. Therefore, we assumed that the HPLM algorithm is appropriate for phenology detection combining Landsat 7 and 8 "overlapping" regions where the temporal resolution could be as fine as 4 days. Land surface phenology in this study was extracted from ETM+ and OLI NDVI ( Figure 2) using the HPLM algorithm. The HPLM divided a vegetation growing cycle into greenup phase and senescence phase using the slope within a moving window of five NDVI observations. In each phase, we fitted NDVI time series using HPLM: where t is time in day of year (DOY); a 1 and a 2 are related to the vegetation growth time; b 1 and b 2 are associated with the rate of plant leaf development; c 1 and c 2 are the amplitude of the NDVI variation; d is the vegetation stress factor; and NDV I b is the background NDVI value. Equation (1) represents the favorable growth condition and Equation (2) denotes vegetation growth suffering from stress. After comparing the fitting goodness of Equations (1) and (2), the better equation was chosen for extracting phenological timing metrics. The phenological transition dates were then identified using the rate of curvature change on the NDVI curve modeled by HPLM [45]. Specifically, the timings corresponding to the extreme values in the rate of curvature change were considered as phenological dates, namely, greenup onset date (GO) and dormancy onset date (DO) (Figure 2).

Evaluation of Phenological Extraction
The confidence level of phenological extraction was measured by the proportion of good-quality observations during a vegetation growing cycle and the fitness of HPLM modeling of the NDVI time series [44]. The proportion of good-quality observations (PGQ) during the growing season was used to qualify if the phenology is appropriately detected. The PGQ can only be calculated after phenological date detections because phenological metrics could be well detected even if missing land surface observations occur during the vegetation dormancy period. In general, high PGQ corresponds to high accuracy in phenology detections. However, it does not necessarily indicate that phenology detections are always incorrect if PGQ is low. The PGQ is quantified as follows: where is the total number of NDVI datapoints during the growing season, and is the number of NDVI datapoints with good quality. PGQ is set to 0 if there are no observations consecutively for more than one month.
The fitness of HPLM modeling was further characterized by comparing modeled values and satellite measurements with cloud-free observations at each pixel. An index of agreement (d) was defined as follows: where denotes the fitness of the model (ranging from 0% to 100%), is the number of NDVI datapoints with good quality, ( ) is the HPLM-modeled NDVI values, ( ) is the observed NDVI values with good quality, and ̅ is the mean value of ( ) . This index is comparable across different vegetated biomes for the assessment of temporal NDVI curve fitting [44,49].
For assessing the reliability of phenological extraction, we calculated the average proportion of good-quality observations (PGQ) and the average index of agreement (d) of HPLM fitting of NDVI curves from 2013 to 2015. Generally speaking, PGQ scores in eastern zones (Zones 1, 2, and 4-8) were lower than in western zones (Zones 9-17) because cloud cover is commonly higher in eastern parts than in western parts of the QTP. The largest PGQ scores appeared in the central area of each zone, where Landsat observed with the highest frequency in a year. In addition, the d values were mostly higher than 80%, indicating that HPLM can effectively fit the NDVI curves ( Figure 3).

Evaluation of Phenological Extraction
The confidence level of phenological extraction was measured by the proportion of good-quality observations during a vegetation growing cycle and the fitness of HPLM modeling of the NDVI time series [44]. The proportion of good-quality observations (PGQ) during the growing season was used to qualify if the phenology is appropriately detected. The PGQ can only be calculated after phenological date detections because phenological metrics could be well detected even if missing land surface observations occur during the vegetation dormancy period. In general, high PGQ corresponds to high accuracy in phenology detections. However, it does not necessarily indicate that phenology detections are always incorrect if PGQ is low. The PGQ is quantified as follows: where T is the total number of NDVI datapoints during the growing season, and n is the number of NDVI datapoints with good quality. PGQ is set to 0 if there are no observations consecutively for more than one month. The fitness of HPLM modeling was further characterized by comparing modeled values and satellite measurements with cloud-free observations at each pixel. An index of agreement (d) was defined as follows: where d denotes the fitness of the model (ranging from 0% to 100%), n is the number of NDVI datapoints with good quality, y (i) is the HPLM-modeled NDVI values, o (i) is the observed NDVI values with good quality, and o is the mean value of o (i) . This index is comparable across different vegetated biomes for the assessment of temporal NDVI curve fitting [44,49]. For assessing the reliability of phenological extraction, we calculated the average proportion of good-quality observations (PGQ) and the average index of agreement (d) of HPLM fitting of NDVI curves from 2013 to 2015. Generally speaking, PGQ scores in eastern zones (Zones 1, 2, and 4-8) were lower than in western zones (Zones 9-17) because cloud cover is commonly higher in eastern parts than in western parts of the QTP. The largest PGQ scores appeared in the central area of each zone, where Landsat observed with the highest frequency in a year. In addition, the d values were mostly higher than 80%, indicating that HPLM can effectively fit the NDVI curves ( Figure 3).

Statistical Analysis of Terrain Effect on LSP
The terrain effects on phenological transition dates were assessed using the high-resolution digital elevation model (DEM) and LSP data. The DEM data at a spatial resolution of 30 m on the QTP were acquired from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global digital elevation model (GDEM) at the USGS EROS Center (United States Geological Survey Earth Resources Observation and Science, http://earthexplorer.usgs.gov/). First, the elevation, slope, and aspect in the 17 Landsat zones were extracted from ASTER DEM using the spatial analysis tools in ArcGIS 10.3 (Figures A2a-c in Appendix A). Aspect was divided into eight classes according to the ArcGIS aspect classification standard in clockwise order, namely, north, northeast, east, southeast, south, southwest, west, and northwest ( Figure A2b). Flat was excluded as it accounted for only a tiny portion of pixels. Three aspects were considered as "shaded", namely, northwest, north, and northeast. These aspects were characterized by less direct sunshine, lower temperature, lower evaporation, and higher moisture [50]. In contrast, the "sunlit" aspects included southeast, south, and southwest. The slopes are larger in the eastern areas (average slope = 19.82° from Zones 1 to 7) than in the western areas (average slope = 10.90° from Zones 8 to 17) ( Figure A2c).
Further, the relationship between the phenological transition metrics and aspect was investigated by removing the impacts from the vegetation type, elevation, and slope. Specifically, the pixels over the 17 zones were classified into different groups according to vegetation types (meadow and steppe), elevations (from 1000 m to 6000 m above mean sea level), and slopes (0° to 70°). Based on one-meter (elevation) and one-degree (slope) intervals, there were up to 5000 × 71 possible combinations (groups) for each vegetation type. In each group, the respective mean phenological dates were computed for all pixels of the eight aspects. The difference between the mean phenological date of a given aspect and the mean phenological date of the eight aspects in a group was defined as the relative phenological date for the given aspect and the given group. By averaging relative phenological dates for a given aspect over all groups, the relative phenological date for a given aspect across the belt transect of the QTP was obtained.
Similarly, the relationship between the phenological transition metrics and slope was analyzed by eliminating the impacts from the vegetation type, elevation, and aspect. Specifically, the pixels with the same vegetation type, elevation, and aspect were classified into one group. Considering elevations from 1000 m to 6000 m and eight classes of aspects, there were up to 5000 × 8 combinations (groups) for a specific vegetation type. In a group, the mean phenological date was computed for all pixels of each slope class, respectively. The difference between the mean phenological date of a given slope class and the mean phenological date of all slope classes in a group was defined as the relative phenological date for the given slope class and the given group. The relative phenological date for a given slope class across the belt transect of the QTP was calculated by averaging relative phenological

Statistical Analysis of Terrain Effect on LSP
The terrain effects on phenological transition dates were assessed using the high-resolution digital elevation model (DEM) and LSP data. The DEM data at a spatial resolution of 30 m on the QTP were acquired from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global digital elevation model (GDEM) at the USGS EROS Center (United States Geological Survey Earth Resources Observation and Science, http://earthexplorer.usgs.gov/). First, the elevation, slope, and aspect in the 17 Landsat zones were extracted from ASTER DEM using the spatial analysis tools in ArcGIS 10.3 ( Figure A2 in Appendix A). Aspect was divided into eight classes according to the ArcGIS aspect classification standard in clockwise order, namely, north, northeast, east, southeast, south, southwest, west, and northwest ( Figure A2b). Flat was excluded as it accounted for only a tiny portion of pixels. Three aspects were considered as "shaded", namely, northwest, north, and northeast. These aspects were characterized by less direct sunshine, lower temperature, lower evaporation, and higher moisture [50]. In contrast, the "sunlit" aspects included southeast, south, and southwest. The slopes are larger in the eastern areas (average slope = 19.82 • from Zones 1 to 7) than in the western areas (average slope = 10.90 • from Zones 8 to 17) ( Figure A2c).
Further, the relationship between the phenological transition metrics and aspect was investigated by removing the impacts from the vegetation type, elevation, and slope. Specifically, the pixels over the 17 zones were classified into different groups according to vegetation types (meadow and steppe), elevations (from 1000 m to 6000 m above mean sea level), and slopes (0 • to 70 • ). Based on one-meter (elevation) and one-degree (slope) intervals, there were up to 5000 × 71 possible combinations (groups) for each vegetation type. In each group, the respective mean phenological dates were computed for all pixels of the eight aspects. The difference between the mean phenological date of a given aspect and the mean phenological date of the eight aspects in a group was defined as the relative phenological date for the given aspect and the given group. By averaging relative phenological dates for a given aspect over all groups, the relative phenological date for a given aspect across the belt transect of the QTP was obtained.
Similarly, the relationship between the phenological transition metrics and slope was analyzed by eliminating the impacts from the vegetation type, elevation, and aspect. Specifically, the pixels with the same vegetation type, elevation, and aspect were classified into one group. Considering elevations from 1000 m to 6000 m and eight classes of aspects, there were up to 5000 × 8 combinations (groups) for a specific vegetation type. In a group, the mean phenological date was computed for all pixels of each slope class, respectively. The difference between the mean phenological date of a given slope class and the mean phenological date of all slope classes in a group was defined as the relative phenological date for the given slope class and the given group. The relative phenological date for a given slope class across the belt transect of the QTP was calculated by averaging relative phenological dates in all possible groups. In order to detect combined effects of slope and aspect on the phenological transition metrics, we also analyzed phenological dates on different slopes of south and north aspects, which represent the typical sunlit and shaded aspects.

Effects of Topographical Parameters on Land Surface Phenology at Local Scales
To reveal detailed spatial patterns of Landsat LSP and its relation to topographical parameters, we extracted elevation, aspect, slope, and greenup onset and dormancy onset dates within a 12 × 12 km steppe area in the 17th Landsat zone. In this area, the elevation varies from 4300 to 6250 m (Figure 4a). The earlier greenup dates (125 to 145 DOY) and later dormancy dates (285 to 320 DOY) (i.e., the longer growing season) appeared in the valley where elevation is less than 4400 m and slope is less than 10 • (Figure 4a,c,e,f). From the piedmont to the snow line near the mountaintop, the greenup onset tended to be delayed at a rate of 0.54 day/100 m (r = 0.64, p < 0.0001), and dormancy onset tended to advance at a rate of −0.46 day/100 m (r = −0.81, p < 0.0001). However, LSP was undetectable in areas with an elevation above 5600 m, where snow generally covered over the entire year (Figure 4d,e,f). Near the snow line, greenup and dormancy onset usually occurred in mid-July (190 to 200 DOY) and mid-September (250 to 260 DOY) with a growing season duration of less than two months. This pattern clearly demonstrates the strong impacts of elevation on vegetation phenology. dates in all possible groups. In order to detect combined effects of slope and aspect on the phenological transition metrics, we also analyzed phenological dates on different slopes of south and north aspects, which represent the typical sunlit and shaded aspects.

Effects of Topographical Parameters on Land Surface Phenology at Local Scales
To reveal detailed spatial patterns of Landsat LSP and its relation to topographical parameters, we extracted elevation, aspect, slope, and greenup onset and dormancy onset dates within a 12 × 12 km steppe area in the 17th Landsat zone. In this area, the elevation varies from 4300 to 6250 m ( Figure  4a). The earlier greenup dates (125 to 145 DOY) and later dormancy dates (285 to 320 DOY) (i.e., the longer growing season) appeared in the valley where elevation is less than 4400 m and slope is less than 10° (Figures 4a,c,   The spatial pattern of aspect also shows the significant impacts on vegetation phenology (Figure 4b,e,f). Taking the Mountains 1, 2, and 3 as examples (Figure 4d), greenup and dormancy onsets on the sunlit slopes (south and southeast) were both about 20-30 days later than those on the shaded slopes (north and northwest). Figure 5 shows spatial patterns of three-year average greenup onset dates and dormancy onset dates over the 17 zones. Greenup onset dates occurred in late March in eastern regions and irregularly shifted westwards with dates as late as July. In contrast, dormancy onset appeared from early August to late November but did not show an obvious gradient change from east to west.  Figure 5 shows spatial patterns of three-year average greenup onset dates and dormancy onset dates over the 17 zones. Greenup onset dates occurred in late March in eastern regions and irregularly shifted westwards with dates as late as July. In contrast, dormancy onset appeared from early August to late November but did not show an obvious gradient change from east to west. With respect to relationship between land surface phenology and elevation, the average greenup onset date showed a delayed tendency (p < 0.0001) from low elevation to high elevation, while the average dormancy onset date indicated an advanced tendency (p < 0.0001) from low elevation to high elevation. The gradients of average greenup onset date and average dormancy onset date along elevation from 1000 m to 6000 m were 1.52 days/100 m and −0.59 days/100 m, respectively. It should be noted that phenological progression with the elevation increase was interrupted in some elevation ranges, particularly around 3400 m. In addition, the dormancy onset dates turned to delay at elevations above 5600 m ( Figure 6). To reveal the relationship between land surface phenology and aspect, we calculated the mean relative phenological date for a given aspect (difference between the mean phenological date on a given aspect and the mean phenological date on all aspects across the belt transect of the QTP) in meadow and steppe areas, respectively. For the timing of greenup onset in the meadow area, the relative greenup onset dates were mostly positive on north and northwest slopes (shaded slopes) during 2013 to 2015, indicating a relatively later greenup onset date. In contrast, the relative greenup onset dates were negative on east and southeast slopes (sunlit slopes), indicating a relatively earlier greenup onset date. Specifically, greenup onset was 1 to 2 days earlier on the partially sunlit slopes With respect to relationship between land surface phenology and elevation, the average greenup onset date showed a delayed tendency (p < 0.0001) from low elevation to high elevation, while the average dormancy onset date indicated an advanced tendency (p < 0.0001) from low elevation to high elevation. The gradients of average greenup onset date and average dormancy onset date along elevation from 1000 m to 6000 m were 1.52 days/100 m and −0.59 days/100 m, respectively. It should be noted that phenological progression with the elevation increase was interrupted in some elevation ranges, particularly around 3400 m. In addition, the dormancy onset dates turned to delay at elevations above 5600 m ( Figure 6).  Figure 5 shows spatial patterns of three-year average greenup onset dates and dormancy onset dates over the 17 zones. Greenup onset dates occurred in late March in eastern regions and irregularly shifted westwards with dates as late as July. In contrast, dormancy onset appeared from early August to late November but did not show an obvious gradient change from east to west. With respect to relationship between land surface phenology and elevation, the average greenup onset date showed a delayed tendency (p < 0.0001) from low elevation to high elevation, while the average dormancy onset date indicated an advanced tendency (p < 0.0001) from low elevation to high elevation. The gradients of average greenup onset date and average dormancy onset date along elevation from 1000 m to 6000 m were 1.52 days/100 m and −0.59 days/100 m, respectively. It should be noted that phenological progression with the elevation increase was interrupted in some elevation ranges, particularly around 3400 m. In addition, the dormancy onset dates turned to delay at elevations above 5600 m ( Figure 6). To reveal the relationship between land surface phenology and aspect, we calculated the mean relative phenological date for a given aspect (difference between the mean phenological date on a given aspect and the mean phenological date on all aspects across the belt transect of the QTP) in meadow and steppe areas, respectively. For the timing of greenup onset in the meadow area, the relative greenup onset dates were mostly positive on north and northwest slopes (shaded slopes) during 2013 to 2015, indicating a relatively later greenup onset date. In contrast, the relative greenup onset dates were negative on east and southeast slopes (sunlit slopes), indicating a relatively earlier greenup onset date. Specifically, greenup onset was 1 to 2 days earlier on the partially sunlit slopes To reveal the relationship between land surface phenology and aspect, we calculated the mean relative phenological date for a given aspect (difference between the mean phenological date on a given aspect and the mean phenological date on all aspects across the belt transect of the QTP) in meadow and steppe areas, respectively. For the timing of greenup onset in the meadow area, the relative greenup onset dates were mostly positive on north and northwest slopes (shaded slopes) during 2013 to 2015, indicating a relatively later greenup onset date. In contrast, the relative greenup onset dates were negative on east and southeast slopes (sunlit slopes), indicating a relatively earlier greenup onset date. Specifically, greenup onset was 1 to 2 days earlier on the partially sunlit slopes than on the shaded slopes (Figure 7a). For the timings of dormancy onset in the meadow area however, the earlier dates occurred on north, northeast, and northwest slopes (shaded slopes), while the later dates appeared on southeast, south, and southwest slopes (sunlit slopes). Overall, dormancy onset was about 2-4 days earlier on the shaded slopes than on the sunlit slopes (Figure 7b). According to a variance analysis and Student's t-test, all of the differences are significant (p < 0.05).

Effects of Topographical Parameters on Land Surface Phenology at Regional Scales
In the steppe area, the positive relative greenup onset dates were mostly found on partially sunlit slopes (southwest and west slopes), while the negative relative greenup onset dates were detected not only on shaded slopes (north and northeast slopes), but also on sunlit slopes (south and southeast slopes) (Figure 7c). Similar to the pattern in the meadow area, dormancy onset on shaded slopes (north, northeast, and northwest) occurred about 3-5 days earlier than that on sunlit slopes (southeast, south, and southwest) (Figure 7d). Based on a variance analysis and Student's t-test, all of the differences are also significant (p < 0.05).
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 17 than on the shaded slopes (Figure 7a). For the timings of dormancy onset in the meadow area however, the earlier dates occurred on north, northeast, and northwest slopes (shaded slopes), while the later dates appeared on southeast, south, and southwest slopes (sunlit slopes). Overall, dormancy onset was about 2-4 days earlier on the shaded slopes than on the sunlit slopes (Figure 7b). According to a variance analysis and Student's t-test, all of the differences are significant (p < 0.05).
In the steppe area, the positive relative greenup onset dates were mostly found on partially sunlit slopes (southwest and west slopes), while the negative relative greenup onset dates were detected not only on shaded slopes (north and northeast slopes), but also on sunlit slopes (south and southeast slopes) (Figure 7c). Similar to the pattern in the meadow area, dormancy onset on shaded slopes (north, northeast, and northwest) occurred about 3-5 days earlier than that on sunlit slopes (southeast, south, and southwest) (Figure 7d). Based on a variance analysis and Student's t-test, all of the differences are also significant (p < 0.05). Considering the effects of the slope on land surface phenology in the meadow area, the mean relative greenup onset date shows nonlinear variation with the slope increase on both north ( Figure  8a) and south (Figure 8c) slopes. The relative greenup date remained relatively stable without a significant trend on slopes less than 30°. However, the relative greenup date correlated significantly positively (p < 0.0001) with slope on slopes more than 30° on both north (Figure 8a) and south ( Figure  8c) slopes. More specifically, mean relative greenup onset date was delayed with the increase of slope (from 0-50°) at rates of 1.15 (on north slope; p < 0.0001) and 0.72 (on south slope; p < 0.0001) days per 10°. In contrast, mean relative dormancy onset date significantly advanced (p < 0.0001) with the increase of slope at rates of 2.74 (on north slope; Figure 8b) and 0.78 (on south slope; Figure 8d) days Considering the effects of the slope on land surface phenology in the meadow area, the mean relative greenup onset date shows nonlinear variation with the slope increase on both north (Figure 8a) and south (Figure 8c) slopes. The relative greenup date remained relatively stable without a significant trend on slopes less than 30 • . However, the relative greenup date correlated significantly positively (p < 0.0001) with slope on slopes more than 30 • on both north (Figure 8a) and south (Figure 8c) slopes. More specifically, mean relative greenup onset date was delayed with the increase of slope (from 0-50 • ) at rates of 1.15 (on north slope; p < 0.0001) and 0.72 (on south slope; p < 0.0001) days per 10 • . In contrast, mean relative dormancy onset date significantly advanced (p < 0.0001) with the increase of slope at rates of 2.74 (on north slope; Figure 8b) and 0.78 (on south slope; Figure 8d) days per 10 • . Overall, as slope increases, greenup and dormancy onset dates change more rapidly on the north slope than on the south slope in the meadow area.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 17 are 3.6 days/10° on the north slope and 5 days/10° on the south slope. The mean relative dormancy onset date is slightly delayed (p > 0.05) with the increase of slope on the north slope (Figure 9b) but significantly delayed (p < 0.0001) with the increase of slope on the south slope (Figure 9d). The mean changing rate on the south slope is 3.44 days/10°. This pattern indicates that greenup and dormancy onset dates changed more rapidly with the increase of slope on the south slope than on the north slope in the steppe area.  Error bars denote the standard deviations of relative phenological dates at given slopes.
In the steppe area, however, the impact pattern of slope was opposite to that in the meadow area. The mean relative greenup onset date significantly shifted to an earlier date with the increase of slope (p < 0.0001) on both north ( Figure 9a) and south (Figure 9c) slopes. The mean changing rates are 3.6 days/10 • on the north slope and 5 days/10 • on the south slope. The mean relative dormancy onset date is slightly delayed (p > 0.05) with the increase of slope on the north slope (Figure 9b) but significantly delayed (p < 0.0001) with the increase of slope on the south slope (Figure 9d). The mean changing rate on the south slope is 3.44 days/10 • . This pattern indicates that greenup and dormancy onset dates changed more rapidly with the increase of slope on the south slope than on the north slope in the steppe area.  north (a,b) and south (c,d) slopes in the meadow area. Dotted lines represent mean relative phenological dates at given slopes (interval = 1°). Error bars denote the standard deviations of relative phenological dates at given slopes.

Discussion
Our research investigated a new approach to detecting vegetation phenology from the Landsat observations in the same year and then quantified the impacts of topographic features at 30 m resolution on the spatial heterogeneity of vegetation phenological metrics. It improved the understanding of vegetation phenology variation on the QTP. Combining Landsat ETM+/OLI observations in overlapping zones of adjacent images offered a high frequency of data for the extraction of land surface phenology. This increase in scene number allowed us to investigate the spatial patterns of vegetation phenology in the same year at a spatial resolution of 30 m in a region spanning 2600 km of the QTP. The reliability of LSP extractions is strongly dependent on the quality of Landsat observations, which is influenced mainly by the number of Landsat overpasses, frequency of cloud contamination, and the Scan Line Corrector (SLC) failure in ETM+. Therefore, to reveal spatial variations of LSP precisely, the evaluation of Landsat observations and phenological extractions is necessary. As the quality of Landsat-derived NDVI data can be measured by PGQ-namely, the higher the PGQ score, the higher the quality of NDVI data-NDVI data in western zones and the central area of each zone should be better than in eastern zones and the edge areas of each zone on the QTP. In this study, we removed NDVI data with the PGQ less than 25% because the LSP with PGQ above 25% is closer to ground-based phenological date observations acquired from the limited ground stations on the QTP. Our study confirmed that PGQ could be used as an important measure of data quality in Landsat-based vegetation phenology extractions. However, PGQ provides an evaluation of overall quality for an EVI2 time series, but it is less likely able to quantify the detection quality of a specific phenological transition date. Indeed, it has been demonstrated that large uncertainties in phenology detections are produced when missing Landsat observations occur around phenological transition dates [48]. Thus, the local proportion of good-quality (LPGQ) observations has been proposed to evaluate the detection quality for each phenological transition date in a vegetation growing cycle [36], which is calculated based on the good-quality observations around the given phenological date. Certainly, LPGQ provides a better evaluation of NDVI data quality around each phenological date in an individual pixel than does PGQ. In this study, however, as topographic impacts on spatial heterogeneity of vegetation phenological metrics were derived by analyzing a large number of pixels at a spatial resolution of 30 m, NDVI data quality around each phenological date is assumed to have only marginal influences on the phenology timing determination.
Based on combined analysis of LSP and elevation, greenup onset is delayed but dormancy onset advances as the elevation increases on the QTP, which is consistent with the results derived from coarse-spatial-resolution satellite data (1-8 km) [10,14,15,19]. However, the changing rates in our study are significantly larger than those in the previous studies (Table 1) because scaling effects may induce larger changing rates at the spatial resolution of 30 m [51]. In other words, the coarser-resolution satellite data may smooth and attenuate the changing rates of LSP. Thus, the elevation dependence of change in land surface phenology detected by Landsat data should be more precise than that derived by AVHRR, MODIS, and SPOT data. Moreover, the elevational gradient of change in greenup onset date was nearly four times as large as that of dormancy onset (Figure 6), which confirms the conclusion that spring phenology variation is very sensitive to air temperature variation [9,52]. By contrast, the smaller elevational gradient of change in dormancy onset date is likely associated with the strong impacts of both temperature and precipitation [15], because precipitation is less dependent on elevation than air temperature. Meanwhile, interruption of phenological progression at about 3400 m might be associated with abrupt elevation change from Zone 1 to Zone 2 ( Figure A2), while the delayed tendency of the dormancy onset date above 5600 m (Figure 6b) may be induced by disturbance of snow cover (see Figure 4d,e,f).
The elevation dependence of vegetation LSP on the QTP has been the research focus for terrain effects; whereas, influences of aspect and slope on LSP have usually been neglected. However, the topographic controls of local vegetation phenology in the mountainous region are integrated by elevation, aspect, and slope. Landsat data with the higher spatial resolution allowed us to detect the impact of terrain on phenology through aspect and slope. To quantify the separate effects of the aspect or slope on the LSP, we used the single variable principle to eliminate the other factors which controlled the LSP. The results confirm that aspect and slope do have strong controls on LSP and the climatic cues are diverse. In the meadow area, greenup onset on shaded slopes was usually later than that on sunlit slopes, which can be explained by lower temperature and later snowmelt timing on the shaded slopes than on the sunlit slopes [50]. Contrarily, the lower temperature on the shaded slopes triggered an earlier occurrence date of dormancy onset, while the higher temperature on the sunlit slopes induced a later occurrence date of dormancy onset (Figure 7a,b). In the steppe area, greenup onset was earlier on both shaded and sunlit slopes (Figure 7c), which indicates that greenup onset may not be mainly affected by air temperature but by precipitation and soil moisture regimes [8]. As the dependence of precipitation on aspect may be smaller than that of temperature in the steppe area where slopes are usually smaller than 30 • (Figures 1 and A2), the differences in greenup onset dates between shaded and sunlit slopes are not significant. However, similar to the dormancy onset response to aspect in the meadow area, the dormancy onset in the steppe area was also earlier on shaded slopes than on sunlit slopes (Figure 7d), which may be induced by lower temperature. Moreover, the effect of aspect on phenological occurrence dates was larger in the steppe area than in the meadow area, especially for dormancy onset. This finding implies that the response of steppe dormancy onset to temperature regimes might be more sensitive than that of meadow dormancy onset on the QTP, possibly due to the lower levels of moisture in the steppe.
With respect to slope effects on LSP, we detected opposite patterns in meadow and steppe areas (Figures 8 and 9). To analyze the climatic attribution of slope effects on LSP, relationships between slope steepness and microclimatic factors should be identified. Unfortunately, there are no specified meteorological observations on different aspects and slopes of the QTP. Thus, we can only roughly evaluate possible climatic cues of slope effects on LSP. According to microclimatic observations on different aspects and slopes of mountainous areas in southern China [53], soil temperature and evaporation decrease with the increase of slope on the north slope but increase on the south slope. Statistical analysis and process-based modeling of long-term phenological observation data on the QTP showed that herbaceous plant green leaves may not grow until the ground has been warmed above 0 • C after snow melt. Hence, snow-melt water is a crucial factor for triggering greenup in addition to temperature in snow-dominated regions of the QTP [9]. That is, more snowfall during late winter and early spring can translate into more soil moisture storage as spring temperature increases, and induce earlier greenup; whereas, less snowfall during late winter and early spring can limit soil moisture availability under rapid spring temperature increase, and force later greenup. In addition, statistical analysis based on remote sensing data on the QTP indicated that the end date of the growing season (EGS) was also controlled by thermal-moisture conditions. Namely, the higher the monthly mean air temperature and the more the monthly precipitation, the later the EGS [15]. Based on the above microclimatic characteristics and plant phenological responses to climatic factors, we attempt to explain slope effects on LSP.
In the meadow area, the range of slope steepness with effective vegetation phenological data was between 0 and 50 degrees. On the south slope, although soil temperature might be high on the steep slope, strong evaporation may induce soil water deficit and limit the advancing effect of high temperature on greenup onset and delaying effect of high temperature on dormancy onset. Due to predominant effects of soil water deficit on plant phenology on the south slope, later greenup onset and earlier dormancy onset occurred on the steep slope (Figure 8c,d). On the north slope, however, lower soil temperature on the steep slope may play a vital role and result in later greenup onset and earlier dormancy onset (Figure 8a,b). In the steppe area, the climate is drier than in the meadow area. The plant community is mainly composed of meso-xerophytes and xerophytes, such as Stipa and Artemisia. As the range of slope steepness with effective vegetation phenological data was less than in the meadow area (between 0 and 40 degrees), the combination of high soil temperature and strong evaporation on the steep south slope may not induce significant soil water shortage. Under the leading effect of the soil temperature gradient change on plant phenology, earlier greenup onset and later dormancy onset appeared on the steep slope (Figure 9c,d). On the north slope, the range of soil temperature and evaporation decrease with the increase of slope might be also less than that in the meadow area, and result in a slightly lower soil temperature but relatively higher soil moisture on the steep slope. Under the leading effect of the soil moisture gradient change on plant phenology, greenup onset shifted to an earlier date and dormancy onset was postponed slightly with the increase of slope (Figure 9a,b).

Conclusions
This study has realized high-precision extraction of vegetation phenology through combining Landsat ETM+/OLI observations in overlapping zones of adjacent images from 2013 to 2015, and revealed terrain effects on land surface phenology along a 2600 km belt transect on the QTP. This approach is of crucial importance for accurately estimating phenological performance under complex terrain conditions. The main conclusions are as follows.
Within a 12 × 12 km steppe area, spatial patterns of greenup onset and dormancy onset dates were controlled mainly by elevation and aspect. At regional scales, the greenup onset date showed a delayed tendency with the increase of elevation at a mean rate of 1.52 days/100 m, while the dormancy onset date indicated an advanced tendency at a mean rate of −0.59 days/100 m. The elevation dependence of change in land surface phenology detected by Landsat data in our study is significantly larger and more precise than those derived from AVHRR, MODIS, and SPOT data in the previous studies. In the meadow area, greenup onset date was usually later but dormancy onset date was commonly earlier on shaded slopes than on sunlit slopes. The main attribution of phenological timing distribution on different aspects might be temperature regimes. By contrast, greenup onset date in the steppe area did not significantly depend on aspects. Precipitation and soil moisture regimes may play a role in causing greenup onset timing distribution on different aspects. However, dormancy onset date in the steppe area indicated a similar response to aspects with that in the meadow area. With regard to the effect of slope on vegetation phenology in the meadow area, greenup onset date was significantly delayed but dormancy onset date significantly advanced with the increase of slope on both north and south slopes. Either delaying rate in greenup onset or advancing rate in dormancy onset was more rapid on the north slope than on the south slope. In the steppe area, however, greenup onset significantly shifted to an earlier date with the increase of slope on both north and south slopes, while dormancy onset was significantly postponed with the increase of slope only on the south slopes. The advancing rate in greenup onset and delaying rate in dormancy onset were more rapid on the south slope than on the north slope. The coordination status of temperature and moisture on different aspects with the increase of slope controls gradient changes of spring and autumn phenology. Funding: This work was funded by the National Natural Science Foundation of China under grant Nos. 41471033 and 41771049, and the scholarship of the China Scholarship Council. surface temperature, and Aster DEM data, and Google Earth for offering the high-resolution images. S.A. also thanks the Geospatial Sciences Center of Excellence at South Dakota State University for providing in-kind support.