Quantitative Assessment of the Impact of Physical and Anthropogenic Factors on Vegetation Spatial-Temporal Variation in Northern Tibet

Abstract

The alpine grassland on the Qinghai-Tibet Plateau covers an area of about 1/3 of China’s total grassland area and plays a crucial role in regulating grassland ecological functions. Both environmental changes and irrational use of the grassland can result in severe grassland degradation in some areas of the Qinghai-Tibet Plateau. However, the magnitude and patterns of the physical and anthropogenic factors in driving grassland variation over northern Tibet remain debatable, and the interactive influences among those factors are still unclear. In this study, we employed a geographical detector model to quantify the primary and interactive impacts of both the physical factors (precipitation, temperature, sunshine duration, soil type, elevation, slope, and aspect) and the anthropogenic factors (population density, road density, residential density, grazing density, per capita GDP, and land use type) on vegetation variation from 2000 to 2015 in northern Tibet. Our results show that the vegetation index in northern Tibet significantly decreased from 2000 to 2015. Overall, the stability of vegetation types was sorted as follows: the alpine scrub > the alpine steppe > the alpine meadow. The physical factors, rather than the anthropogenic factors, have been the primary driving factors for vegetation dynamics in northern Tibet. Specifically, meteorological factors best explained the alpine meadow and alpine steppe variation. Precipitation was the key factor that influenced the alpine meadow variation, whereas temperature was the key factor that contributed to the alpine steppe variation. The anthropogenic factors, such as population density, grazing density and per capita GDP, influenced the alpine scrub variation most. The influence of population density is highly similar to that of grazing density, which may provide convenient access to simplify the study of the anthropogenic activities in the Tibet plateau. The interactions between the driving factors had larger effects on vegetation than any single factor. In the alpine meadow, the interaction between precipitation and temperature can explain 44.6% of the vegetation variation. In the alpine scrub, the interaction between temperature and GDP was the highest, accounting for 27.5% of vegetation variation. For the alpine steppe, the interaction between soil type and population density can explain 29.4% of the vegetation variation. The highest value of vegetation degradation occurred in the range of 448–469 mm rainfall in the alpine meadow, 0.61–1.23 people/km² in the alpine scrub and –0.83–0.15 °C in the alpine steppe, respectively. These findings could contribute to a better understanding of degradation prevention and sustainable development of the alpine grassland ecosystem in northern Tibet.

1. Introduction

As the main terrestrial ecosystem in the Qinghai-Tibet Plateau, grassland ecosystems are vulnerable to external disturbances and sensitive to environmental changes [1]. In recent years, grassland ecosystems in the Qinghai-Tibet Plateau have faced different degrees of degradation [2]. Both physical and anthropogenic factors (e.g., climate, soil type, topography, population, grazing density and land use [3,4,5,6,7,8]) are considered major driving forces of vegetation variation. Previous studies have analyzed the individual contribution of each physical or anthropogenic factor to variation in vegetation cover through the analysis of correlation and regression coefficients [9,10]. However, the primary and potential interactive or combined impacts of both physical and anthropogenic factors on vegetation variation still remain unclear.

A large number of studies have dedicated enormous efforts to analyze the causes of vegetation variation. In particular, numerous studies have analyzed the relationships between vegetation index and meteorological factors [11,12,13,14,15]. For example, precipitation and temperature have been reported as the main factors affecting variations in grassland vegetation cover in the Qinghai-Tibet plateau [16,17,18]. Some researches have denoted that environmental humidity conditions (precipitation, vapor pressure, and relative humidity) dominate the vegetation index variation over the Tibetan plateau [19]. Still, other studies have indicated that physical factors, such as soil types, solar radiation and elevation, play important roles in vegetation cover change. Soil types determine different water retention and venting properties, and are found to significantly modulate the vegetation trends in dry land [20]. Solar radiation is the dominant controller of plant photosynthesis, and therefore has a huge impact on vegetation productivity in the Qinghai-Tibet plateau [21,22]. With the regard to elevation, some research has shown that vegetation is vulnerable at low altitudes due to the severe influence of human activities [23], while some has found that high-altitude areas in the Tibetan Plateau are especially vulnerable to the degradation process, due to the increased risk of wind and rain erosion [24]. Compared with above-mentioned studies on physical factors, studies on anthropogenic factors are much fewer. Some of the anthropogenic studies have shown that grazing activities lead to severe grassland degradation [5,6,8]. The closer an area is to roads or highways, the more susceptible its vegetation is [25,26]. Previous studies have investigated the individual effect of physical factors or anthropogenic factors on vegetation change. However, the magnitude and patterns of physical and anthropogenic factors in driving grassland variation still remain debatable when these environmental factors are examined together. Most importantly, the combined effects of physical and anthropogenic factors on vegetation change have not been assessed integrally and interactively.

Various models can be used to identify and quantify the physical and anthropogenic factors of vegetation variation. The RESTREND method, proposed by Evans and Geerken [27], is mainly used to distinguish the influence of climate and human activities on vegetation cover change in arid and semi-arid regions. Some mathematical models, such as correlation analysis and multiple regression analysis, can be used to measure the influence of factors on vegetation variation through statistical parameters [10,28]. In addition, the spatial analysis functions of a Geographical Information System (GIS), such as buffer analysis, overlay analysis, and zonal stats, are used to assess the effects of driving factors on vegetation variation [29,30]. However, deficiencies still exist in these current and commonly used research methods, as mentioned above. For example, the RESTREND method is generally appropriate for regions where precipitation is the main limiting factor of vegetation growth. Moreover, some categorical variables, playing important roles in vegetation growth, such as soil type and aspect, cannot be quantified with these methods, and the interaction among influencing factors cannot be measured either.

The geographical detector model, a statistical method proposed by Wang et al, is used to measure both individual and interactive impacts of driving factors on environmental risks based on spatial variation analysis [31,32]. This model has made it possible to reveal the driving factors behind spatial stratified heterogeneity and to quantify the interaction between factors without a linear hypothesis [33,34]. It has been applied in different research fields, such as landscape ecological patterns and health risk assessment [33,35]. With this application, the geographical detector model can be further used as an effective method to explore the spatial correlations between environment factors and vegetation variables.

Northern Tibet, located in the hinterland of the Qinghai-Tibet Plateau, is an important worldwide grassland ecosystem [22,36]. The particularity in its geographical environment and climatic conditions has made northern Tibet a typical region with an extremely vulnerable ecological environment. In this study, we used linear regression analysis and the geographical detector model to quantify the driving factors of vegetation variation from 2000 to 2015 in northern Tibet. We selected both physical factors (precipitation, temperature, sunshine duration, soil type, elevation, slope, and aspect) and anthropogenic factors (population density, road density, residential density, grazing density, per capita GDP, and land use type) for analysis. Our goals were divided into four parts: (1) to identify spatial and temporal trends in vegetation variation in northern Tibet from 2000 to 2015; (2) to quantify the factors, both physical and anthropogenic, which may affect vegetation variation; (3) to reveal how these factors interact to influence vegetation variation; (4) and to find out the areas that are sensitive to vegetation variation.

2. Materials and Methods

2.1. Study Area

Northern Tibet is located in the north part of the Tibetan Autonomous Region (29°55’–36°30’ N, 83°55’–95°05’ E). Northern Tibet adjoins Kunlun Mountains and Danggula Mountains to the northwest, the “Three-River Headwaters” region to the east, and Ali Plateau to the west. It covers approximately 369,000 km², with a mean elevation of over 4500 m. Northern Tibet belongs to the typical sub-frigid zone, where the primary climatic features are dryness and coldness. The annual mean temperature ranges between –28.6 to 9.4 °C, and the annual precipitation varies from 100 mm in the northwest to 800 mm in the southeast during 2000–2015 (interpolated through meteorological station data from the China Meteorological Data Service Center (CMDC), (Figure 2). The growing season typically extends from May to September [15]. About 80% of the total precipitation occurs in this period because of atmospheric circulation and local terrain. Natural vegetation types mainly include alpine shrub, alpine steppe, and alpine meadow (Figure 1). Northern Tibet consists of 11 counties: Nagqu, Nyainrong, Amdo, Baingoin, Xainza, Nyima, Lhari, Biru, Baqen, Sog, and Shuanghu. Except for Biru and Sog, which are semi-agricultural and semi-pasture counties, the nine remaining counties are pasture areas. Animal husbandry on the grassland is the mainstay of the region’s economy. Alpine grassland, the most significant ecosystem in this area, occupies about 90% of the total area in northern Tibet and provides important subsistence for local people [37].

2.2. Datasets

2.2.1. NDVI Time Series

The NDVI data were extracted from the terra Moderate Resolution Imaging Spectro-radiometer (MODIS) Vegetation Indices (MOD13Q1) Version 6 data, with a spatial resolution of 250 m and a temporal resolution of 16 days [38]. Four MODIS tiles, namely h24v05, h25v05, h25v06, and h26v05, were required for the study area. These MODIS tiles were mosaicked and then projected to Albers Conical Equal Area projection through MODIS Reprojection Tool (MRT) software. Then, the monthly NDVI value was calculated with the Maximum Value Composites (MVC) method, which further eliminated errors caused by cloud, atmosphere, sun angle, as well as other interference [39]. The monthly NDVI data from May to September were averaged to yield the growing season NDVI. Finally, to exclude the effects of non-vegetation factors, we extracted areas with average NDVI greater than 0.1 for further analysis [40].

2.2.2. Physical Factors Influencing Vegetation Variation

Meteorological data, including monthly temperature, precipitation and sunshine duration datasets during the period of 2000–2015, were collected from the China Meteorological Data Service Center (CMDC). In total, original data from 39 meteorological stations were interpolated to continuous surface data with a 0.01-degree spatial resolution by using the thin-plate smoothing spline algorithm implemented in ANUSPLIN interpolation software [41].

Topographic data, including elevation, slope, and aspect, were extracted from a digital model (DEM). The DEM was obtained from Shuttle Radar Topography Mission (SRTM) terrain product data with a spatial resolution of 90 m.

Soil type data were extracted from the Soil Map of China (1:1,000,000). The dataset was based on the traditional ‘soil generation classification (CSGC)’, and provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn).

Vegetation type data were obtained from the Vegetation Map of China (1:1,000,000), and provided by the Data Center for Resources and Environmental Sciences, and the Chinese Academy of Sciences (RESDC) (http://www.resdc.cn).

2.2.3. Anthropogenic Factors Influencing Vegetation Variation

Population data, gross domestic product data, and livestock data were acquired from the official statistical surveys published in the Xizang Statistic Yearbook with a county resolution. Population density was calculated by dividing the population of each county by the area of the county. The per capita GDP was obtained by dividing the gross domestic product of the county by the corresponding population. To quantify the grazing density within each county, animal number was standardized to sheep units. Generally, a large animal, such as a horse or a yak, is equivalent to 5 sheep units, and a sheep or a goat equals 1 sheep unit [42].

Spatial Distribution Map of Land Use Types (1:1,000,000) was obtained from the land use status database in China (http://www.resdc.cn). Data production was based on Landsat TM/ETM remote sensing images as the main data source, and then generated by manual visual interpretation.

Road and resident distribution data of the study area were extracted from Open Street Map (OSM) data. With the aid of Kernel Density tool of ArcGIS software, the spatial distribution density data of road and residential areas were obtained.

2.3. Analysis

2.3.1. Mann Kendall (MK) Test

The non-parametric Mann-Kendall, proposed by Mann and Kendall [43,44], is used to detect statistically significant trends in NDVI and continuous environmental variables time series. It does not require the sample to follow a certain distribution, and the test results are basically not interfered by a few abnormal values [45]. The Mann-Kendall test statistic Z is given by:

$$\mathrm{Z}=\{\begin{array}{c}\frac{S-1}{\sqrt{Var\left(S\right)}},\hspace{1em}if S>0\\ 0,\hspace{1em}if S=0\\ \frac{S+1}{\sqrt{Var\left(S\right)}},\hspace{1em}if S<0\end{array}$$

(1)

$$\mathrm{S}={\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{j=i+1}^{n}sgn\left(x_j-x_i\right)$$

(2)

$$sgn\left(x_j-x_i\right)=\{\begin{array}{c}1,\hspace{1em}if x_j-x_i>0\\ 0,\hspace{1em}if x_j-x_i=0\\ -1,\hspace{1em}if x_j-x_i<0\end{array}$$

(3)

$$\mathrm{Var}\left(\mathrm{S}\right)=\frac{n\left(n-1\right)\left(2n+5\right)-{{\displaystyle \sum }}_{i=1}^m{t}_i\left(t_i-1\right)\left(2t_i+5\right)}{18}$$

(4)

where n is the length of time series, m is the number of tied groups, $x_i$, $x_j$ and $t_i$ are the values in time series I, j (j > i) and the number of tied groups in the time series, respectively. A tied group is a set of sample data which have the same value.

The positive values of Z indicate increasing trends, while negative values show decreasing trends. The significance test is estimated by comparing Z and $Z_{\alpha }$. When –$Z_{\alpha }$ < Z < $Z_{\alpha }$, the statistic Z is not statistically significant. When Z > $Z_{\alpha }$ or Z < −$Z_{\alpha }$, Z is statistically significant. In this study, significant level $\alpha$ = 0.05 was used. At the 5% significant level, $Z_{\alpha }$ equals 1.96 correspondingly.

The sequential Mann-Kendall test is used to statistically assess the shift in trends at a regional scale and determine the beginning of a significant trend in a time series [43]. There are two series derived from the test, a progressive one $U{F}_K$ and a backward one, −𝐾. The calculation process of the Mann-Kendall test is as follows:

Firstly, Setting the time series as $y_1$, $y_2$, ... $y_n$, $m_i$ is the cumulative number when sample $y_i$ is greater than sample $y_j$ ($1\le j\le i$). Define statistics:

$$d_k={\displaystyle \sum }_{i=1}^k{m}_i,\left(2\le k\le \mathrm{n}\right)$$

(5)

Then, assuming that the time series are random and independent, the mean and variance of $d_k$ are calculated as:

$$\mathrm{E}\left(d_k\right)=\frac{k\left(k-1\right)}{4}$$

(6)

$$\mathrm{var}\left(d_k\right)=\frac{k\left(k-1\right)\left(2k+5\right)}{72}$$

(7)

At last, a sequential forward series $U{F}_K$ value can be calculated via standardizing $d_k$:

$$U{F}_K=\frac{d_k-E\left(d_k\right)}{\sqrt{var\left(d_k\right)}}$$

(8)

Similarly, this method is applied to calculate backward series $U{B}_K$. If the value of $U{F}_K$ is greater than 0, the sequence shows an upward trend. If it is less than zero, it indicates a downward trend. When the value of $U{F}_K$ exceeds the confidence interval of $\pm 1.96$, it indicates a statistically significant upward or downward trend at the 95% significant level. If the intersection points of $U{F}_K$ and $U{B}_K$ curves appear within the confidence interval of $\pm 1.96$, it is the statistically significant trend shift point at the 95% significant level.

2.3.2. Trend Analysis

We used unary linear regression to detect the change rate of vegetation index NDVI and environment factors at each pixel, and the change rate is further input into the geographical detector model as the dependent variable. The change rate for time series x at pixel scale was calculated as

$$\mathrm{S}=\frac{n\times {{\displaystyle \sum }}_{i=1}^{n}i\times M_{x,i}-{{\displaystyle \sum }}_{i=1}^{n}i\times {{\displaystyle \sum }}_{i=1}^n{M}_{x,i}}{n\times {{\displaystyle \sum }}_{i=1}^n{i}^2-{\left({{\displaystyle \sum }}_{i=1}^{n}i\right)}^2}$$

(9)

where, S is the slope of the time series x; i refers to the ordinal number 1, 2, 3, …, n corresponding to the year from 2000 to 2015; M_x,i is the mean x in the i-th year; when S > 0, the time series x shows an increasing trend, when S < 0, time series x shows a decreasing trend.

2.3.3. Geographical Detector Model

As a method to test spatial heterogeneity and to quantify the impacts of driving factors, the geographical detector model has been widely used in landscape ecological pattern, health risk assessment and the relationships between human activities and geographical conditions [35,46,47,48,49]. We assumed that the dependent variable would exhibit a similar spatial distribution, as the independent variable, if the independent variable contributes to the dependent variable [31]. Specifically, we used this model to test the consistence of spatial patterns between dependent variable y (NDVI tendency, slope of NDVI trend) and independent variable x (the factors influencing the variation of NDVI tendency), and therefore to find the essential driving forces.

The geographical detectors consist of four parts, namely factor detector, interactive detector, risk detector, and ecological detector.

The factor detector is used to quantify the contribution of impact factors. The q value implies the degree to which factor x explains the heterogeneity of y.

$$q_H=1-\frac{1}{N*{\partial }^2}{\displaystyle \sum }_{i=1}^L{N}_i*{\partial }_i^2$$

(10)

where, H represents variable y; N represents the number of total samples; and ${\partial }^2$ represents the variance of H over the whole area. The study area is stratified into L stratums, denoted by i = 1, 2, ……, L, based on the optimal discretization methods listed below for each continuous impactor. $N_i$ represents the number in the i-th sub-region of one impact factor stratum; and ${\partial }_i^2$ represents the variance of H over the i-th sub-region of the stratum. The range of q value is 0 to 1. The larger the value of q is, the greater impact factor x makes on the heterogeneity of y. If factor x completely explains the spatial heterogeneity of y, q value equals 1; if it is totally unrelated to the spatial heterogeneity of y, q value equals 0.

The ecological detector is used to compare whether the impact of two influencing factors on the spatial distribution of vegetation tendency is significantly different or not, using the F test to measure [32]. The interactive detector, which is an advance in geographical detectors over other statistical methods, compares the q values of single factors and of interactive factors, and explores whether the two factors, when taken together, weaken or enhance one another, or have completely independent effects on a research subject. The judgment method of the interaction mode is shown in

Table 1. Types of interaction between two covariates.

Interaction	Judgment Criteria
Enhance	q(X1∩X2) > q(X1) or q(X2)
Enhance, bivariate	q(X1∩X2) > q(X1) and q(X2)
Enhance, nonlinear	q(X1∩X2) > q(X1) + q(X2)
Weaken	q(X1∩X2) < q(X1) + q(X2)
Weaken, univariate	q(X1∩X2) < q(X1) or q(X2)
Weaken, nonlinear	q(X1∩X2) < q(X1) and q(X2)
Independent	q(X1∩X2) = q(X1) + q(X2)

X1 and X2 represent the influence factors of vegetation degradation. The symbol “∩” denotes the interaction between X1 and X2.

. The risk detector is used to identify vulnerable areas under potential risks in determinant strata.

2.3.4. Discretization Method

The geographical detector model is adapted to discrete factors, and can analyze continuous impact factors through appropriate discretization methods to transform the continuous data into discrete data. In this study, several impact factors, such as post-interpolation meteorological factors, elevation, slope and anthropogenic factors, are continuous variables. To avoid randomness and subjectivity of user-defined discretization, we selected the optimal discretization method from five unsupervised discretization methods, namely Equal Interval (EI), Natural Break (NB), Quantile (QU), Geometrical Interval (GI) and Standard Deviation (SD) [50]. q values are used to assess the effectiveness of different discretization methods, for which can reflect the relationships between NDVI tendency and potential driving factors, as described in 2.3.2. Here, temperature is taken as an example to explain the discretization process. First, assigning numbers of cut points. We set the number of discretion intervals as 2–8, and complete the preliminary classification of temperature based on 5 methods and 7 levels of each method. Then, extracting multi-values of the NDVI tendency layer and temperature layers of 35 discretization results. Finally, calculating and comparing the q values of temperature in each classification case. The optimal discretization method and number of intervals are determined when the q value is relatively high. The same discretization process applies to the other continuous variables. The q values of all continuous impact factors in each preliminary classification cases are shown in Appendix A. On a regional scale, the optimal discretization results of physical and anthropogenic factors’ spatial distributions are shown in Figure 2.

Statistical analyses, such as linear model, quadratic model, growth model, logarithm model, cubic model, and exponential model, were conducted to explore the linear and nonlinear relationships between main influencing factors and NDVI tendency by PASW Statistics 18.0. Regression coefficients ($b_0$, $b_1$, $b_2$, and $b_3$) and the goodness of fit (R²) were calculated and compared to determine which model to use.

$$\mathrm{y}=b_0+b_1*x+b_2*x^2+b_3*x^3$$

(11)

$$\mathrm{y}=e^{b_0+b_1*x}$$

(12)

$$\mathrm{y}=b_0+b_1*In\left(x\right)$$

(13)

$$\mathrm{y}=b_0*e^{b_1*x}$$

(14)

$$R^2=\frac{SSR}{SST}=1-\frac{SSE}{SST}=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(15)

where, $b_0$, $b_1$, $b_2$ and $b_3$ are regression coefficients. y represents values of NDVI, x represents values of main environment factors. Each $y_i$ associated with a predicted (or modeled) value $\widehat{y_i}$. $\overline{y}$ is the mean value of the observed data y. SST represents the sum of the squares of the difference of the dependent variable and its mean value; SSR represents the sum of squares of regressions, which is the sum of squared deviations that reflect the degree of correlation between independent and dependent variables; SSE represents sum of squared residuals, which is a measure of the discrepancy between the data and an estimation model.

3. Results

3.1. Trend and Trend Shift Analysis of NDVI and Environment Factors in Northern Tibet

From the southeast to the northwest of northern Tibet, the 16-year average NDVI had been decreasing in growing season (Figure 3). This was consistent with the regional distribution of temperature and precipitation. The northwestern counties (Nyima and Shuanghu counties) had the lowest mean NDVI value (0.16) in the whole region. For the central counties (Amdo, Shenza, and Baingoin counties) the average annual NDVI value was 0.25. In the southeastern counties (Nagqu, Nyainrong, Lhari, Biru, Sog, and Baqing counties), with an average value of 0.45 and maximum value of 0.52, the mean NDVI was higher than those of the other regions.

The spatial distribution of the slope of vegetation and environment factors trend estimated by ULR (Figure 4a–d) and corresponding significant level (P < 0.05) estimated by the MK method (Figure 4e–h) are basically consistent. Areas with recovery tendency were mainly distributed in the north and southeast part of northern Tibet, accounting for 46% (ULR method) and 49% (MK method) of alpine grassland in this region (Figure 4a,e). The southwest and central regions, which account for 54% (ULR method) and 51% (MK method) of the whole region, have already undergone a degradation tendency. In most parts of northern Tibet, the inter-annual NDVI of the alpine grassland experienced small fluctuations. Specifically, significant recovery trends mainly occurred in the north part of northern Tibet, accounting for 9% (ULR method) and 13% (MK method) of the whole area, whereas 10% (ULR method) and 11% (MK method) of the grassland showed significant degradation trends, mainly occurring in the southern part of Northern Tibet. Nagqu county (the capital of Northern Tibet, located in the south of the study area) and its surrounding areas underwent degradation or even significant degradation. Precipitation showed a decreasing trend in most parts of the study area, and declined most dramatically and significantly in the northwest (Figure 4b,f). Areas with increasing trend of precipitation were only distributed in the east part. Temperature showed an increasing trend nearly in the whole region during 2000–2015. The regions with a significant increasing trend (P < 0.05) were mainly distributed in the northeast of northern Tibet (Figure 4c,g). In most of northern Tibet, sunshine duration exhibited a significant decreasing trend (P < 0.05) (Figure 4d,h).

From 2000 to 2015, there have been decreasing trends in the average annual NDVI in the whole region and in the alpine meadow since 2004 (Figure 5a). A decreasing trend shift was found in 2005 and the trend became significant around 2007 in the whole region and in the alpine meadow (P < 0.05). Precipitation had exhibited a decreasing trend since 2006 and reached a significant level (P < 0.05) around 2013 in the whole region, in the alpine scrub and in the alpine steppe (Figure 5b). In particular, there was a decreasing trend shift around 2012 in the alpine scrub. Temperature had been increasing during the whole study period, and the increasing trend was not significant except in the period of 2005–2009 (Figure 5c). We can see that sunshine duration had exhibited a decreasing trend since 2006 in all vegetation types (Figure 5d). The $U{F}_K$ and $U{B}_K$ curves for sunshine duration intersected at 2006, and the intersected point was statistically significant, indicating that a decreasing trend shift began in 2006 and became significant around 2011. Population density across northern Tibet had exhibited an increasing trend since 2002 and became statistically significant around 2005 (P < 0.05) (Figure 5e). Grazing density in northern Tibet had increased first and then decreased since 2012 (Figure 5f).

3.2. Assessing the Impact of Individual Factors on Vegetation Variation

The effects of physical and anthropogenic factors on NDVI tendency differed in three types of grassland. Moreover, the effect of the same individual factor on NDVI of various vegetation types was different (Figure 6). The q values for the alpine meadow were relatively high, which indicated that NDVI distribution was more sensitive in the alpine meadow, than in the other two types of grasslands (

Table 2. Main influencing factors and corresponding q values in different vegetation types. Factor and ecological detector models were used to determine the most powerful factors affecting the spatial distribution of NDVI variation.

Grassland Type	Influencing Factors
Meadow	Precipitation	Grazing density	Population density	Temperature	Per capita GDP	Slope
q	0.282 **	0.205 **	0.201 **	0.165 **	0.157 **	0.141 **
Scrub	Population density	Grazing density	Per capita GDP	Precipitation	Residential density	Temperature
q	0.101 **	0.096 **	0.087 **	0.076 **	0.074 **	0.058 **
Steppe	Temperature	Population density	Per capita GDP	Grazing density	Slope	Road density
q	0.118 **	0.115 **	0.112 **	0.101 **	0.072 **	0.055 **

q values denote the determination power of driving factors on vegetation variation. ** P < 0.01, * P < 0.05.

). The main factor that determined vegetation degradation in the alpine meadow was precipitation (0.282); in the alpine scrub, the main factor was population density (0.101); in the alpine steppe, it was temperature (0.118). Anthropogenic factors, such as population density, grazing density and per capita GDP, influenced all vegetation types hugely. q values and order of them were very close. In the alpine meadow, grazing density (0.205) and population density (0.201) ranked second and third, respectively. In the alpine scrub areas, population density (0.101), grazing density (0.096) and per capita GDP (0.087) ranked top three. In the alpine steppe, population density (0.115), per capita GDP (0.112) and grazing density (0.101) ranked second, third and fourth, respectively.

Most of the physical factors and the tendency of NDVI had better simulation results, with larger r square, compared with anthropogenic factors in the three vegetation types (Figure 7). The most serious vegetation degradation appeared in the areas with fewer rainfalls in the alpine meadow, showing that the tendency of NDVI decreased first and then increased on the spatial scale with the increase of precipitation. In terms of temperature, the most serious vegetation degradation appeared in the areas with higher temperature for the whole region, showing that the tendency of NDVI decreased first and then increased with the increase of temperature in regions of the alpine meadow and the alpine steppe. By contrast, NDVI tendency was significantly negatively correlated with temperature in the alpine scrub (slope = –0.00013, P < 0.01). It is worth noting that terrain slope displayed a significantly positive correlation with NDVI tendency in the alpine meadow (slope = 0.0011, P < 0.05), while in the alpine steppe, terrain slope showed a marginally significantly negative correlation (slope = –0.0001, P = 0.06). In terms of anthropogenic factors, residential density in the alpine scrub, per capita GDP and road density in the alpine steppe, it showed that, generally, the higher each density was, the more serious the vegetation degradation was. However, the alpine scrub presented the opposite result: the smaller the population density was, the more serious the vegetation degradation was.

3.3. Interaction between Factors That Influence Vegetation Variation

All interactive q values, divided into two categories: double factor enhancement and nonlinear enhancement, were greater than the highest q value of any of their individual factors. This indicates that the interaction between factors played an essential role in vegetation variation. In the alpine meadow, the interaction of precipitation with other physical factors was the most significant (

Table 3. The dominant interactions between two impact factors in different vegetation types. The interactive detector was used to check whether two impact factors work independently or interactively.

	Alpine Meadow	Alpine Scrub	Alpine Steppe
1st	Precipitation ∩ Temperature	Temperature ∩ GDP	Soil ∩ Population
q	0.446	0.275	0.294
2nd	Precipitation ∩ Elevation	Sunshine ∩ GDP	Soil ∩ GDP
q	0.443	0.271	0.289
3rd	Precipitation ∩ Sunshine	Sunshine ∩ Population	Soil ∩ Temperature
q	0.416	0.262	0.288

q values denote the determination power of driving factors on vegetation variation.

). Particularly, the interaction between precipitation and temperature can explain 44.6% of the vegetation variation in the alpine meadow. In the alpine scrub, the interaction between temperature and GDP was the highest, which accounted for a 27.5% change of vegetation variation. For the alpine steppe, the top three interactions were the soil type associated with other influencing factors. For example, the interaction between soil and population can explain a 29.4% change of the vegetation variation. In these three vegetation types, the interaction effects of the main factors were greater in the alpine meadow than those in the other two vegetation types.

3.4. Identification of Areas Vulnerable to Vegetation Degradation

The spatial distribution of each vegetation type that is vulnerable to vegetation degradation was significantly different (Figure 8). In the alpine meadow, the risk of vegetation degradation reached the maximum in areas where annual average precipitation levels were between 448 to 469 mm. However, in the alpine scrub, the highest risk of vegetation degradation was observed in the regions with annual average population density between 0.61–1.23 people/km², and in the alpine steppe, it was the areas with annual average temperatures between –0.83 to –0.15 °C.

4. Discussion

The results show that NDVI decreased from 2000 to 2015 across northern Tibet. Previous studies have shown that NDVI tendency varied with the different spatio-temporal scale on the Tibet plateau [10,51]. The inconsistency of NDVI change rate may result from the spatially heterogeneous of climate, land use policies as well as the uncertainty of remote sensing data. The state of climate factors, such as temperature, precipitation and solar radiation varied in ecosystem zones, where the vegetation composition and humidity differed. The state has implemented a pilot project of “returning grazing to grassland” in Tibet since 2004, whose layout and implementation differed on a regional scale. For example, research has shown that the response degree to this “returning grazing to grassland” project in the alpine meadow was obviously higher than that in the alpine grassland [52]. Additionally, both the scale effect and the uncertainty influence the accuracy of the data, causing the inconsistency results of NDVI tendency in studies. Vegetation parameters, such as NDVI, leaf area index (LAI) [53,54], are needed for validation in future studies. The multi-annual average NDVI from our studies shows the same spatial distribution as previous research with a declining trend from the southeast to the northwest of northern Tibet [55], which is consistent with the water and heat distribution pattern in geographic space. Nagqu county and its surrounding areas are well supplied with a good condition of water resource and temperature. However, these areas have been undergoing significant degradation due to the fast society and economy development. Different grassland types’ responses to environment factors differed. NDVI variation in the alpine meadow was generally more sensitive to influencing factors than that in the alpine steppe. This result was similar with the previous studies [16,56]. The alpine steppe in the Tibet plateau is cold-resistant and drought-tolerant, where perennial dense clumps, rhizomes, and shrubs are dominating vegetation. By contrast, the alpine meadow has a relatively simple plant community composition dominated only by perennial herbaceous species. Also, the alpine meadow has high population density; therefore, it is more vulnerable to environmental changes [57].

The dominant driving factors of vegetation variation differed in each vegetation type of northern Tibet, primarily because the areas with diverse vegetation types can present different environmental conditions, and also because the vegetation with specific physiological characteristics has different interactions with environmental factors. Specifically, precipitation was the key factor that influenced the alpine meadow variation, whereas temperature was the key factor that contributed to the alpine steppe variation. Previous studies have yielded different results on the factors’ effects on vegetation variation. For example, Sun and Ding suggested that mean annual precipitation had the greatest effect on the alpine meadow [16,56,58]. Shen et al. found that no significant correlation had been found between the maximum vegetation index and temperature or precipitation in the growing season across the Tibetan Plateau [19]. Obviously, there are still some controversies about the driving factors of different vegetation types. Therefore, further research is needed to explore the relationships between meteorological factors and different vegetation types. Statistical analysis showed that the areas under higher temperature and lower precipitation had higher risk of vegetation degradation. This result agrees well with the findings of Zhang and Cui [55,59]. Increased temperatures and reduced precipitation will worsen atmospheric drought and decrease relative humidity. Hence, this causes a climate drought, which is not conducive to the improvement and maintenance of the grassland ecosystem productivity in the region. The statistical analysis between NDVI tendency and terrain slope indicates that the flatter the terrain in the alpine meadow is, and the steeper the area in the alpine steppe is—the greater the level of vegetation degradation is. These two opposite results may be attributed to the effects of terrain slope on water use efficiency. Where the slope of the alpine steppe is steep, the soil moisture content is low and is not conducive to the growth of vegetation, and the vegetation degradation there is more serious. For the alpine meadow, water is not a limiting factor for plant growth [57]. Instead, human activities are intensive in flat areas. Therefore, vegetation degradation is more serious there.

The anthropogenic factors, such as population density, grazing density and per capita GDP, influenced the alpine scrub variation greatly. Particularly, these factors influenced the alpine meadow and the alpine steppe hugely. The higher the residential density in the alpine scrub is, or the higher per capita GDP or road density in the alpine steppe is, the more serious the vegetation degradation is. This result is parallel with a previous study showing that the degradation of vegetation areas closer to the highway was more severe [30]. We also found that the q values and the order of population density, grazing density and per capita GDP are very close, although the method we used excluded the influence of multiple collinearity among the factors. It indicated that the relationships between the anthropogenic factors are closer than those between the physical factors. Actually, the interference of the anthropogenic factors in the Tibet plateau is not as complicated as that in other livable and developed places. Grazing is the primary factor affecting the ecological environment on the plateau [5,6,8]. Where there are a large number of livestock units, there is often large population density and high social productivity. Therefore, the synthetic influence index of the anthropogenic factors can be considered to simplify the study of quantifying effects of the anthropogenic activities.

The complexity of the geographical processes suggests that multiple factors of environment interactively affect vegetation degradation or recovery. The results of the interactive detector show that the interaction between factors can greatly enhance the explanatory power of vegetation variation. Precipitation was the most significant factor that affected the alpine meadow variation, and it also played an important role in the interaction of factors. These results highlight the importance of precipitation in the alpine meadow. Our results denote that precipitation shows a decreasing trend, temperature shows an increasing trend, and sunshine duration shows a significant decreasing trend on both the temporal and spatial scales (Figure 4b–d), (Figure 5b–d). Previous research has shown that precipitation is positively correlated with vegetation growth on a global scale [60]. The interaction between temperature and precipitation enlarges the impact on vegetation variation because of the decrease of precipitation and the acceleration of evapotranspiration by the increase of temperature, which results in water stress. Altitude is directly related to temperature, precipitation, soil organic matter and vegetation coverage [61,62], and when the effect of precipitation is superimposed with altitude, it has a greater impact on vegetation. Sunshine duration can be used to estimate photosynthetically active radiation (PAR) [63]. PAR refers to the solar spectral components that can be photosynthesized by green plants. It mainly affects the growth and development of plants by controlling photosynthesis [64]. Therefore, the decrease of sunshine duration will aggravate the adverse effects on vegetation caused by the decrease of precipitation. The explanatory power of soil on the alpine steppe variation was greatly enhanced by the interaction with the anthropogenic factors. The alpine steppe has small vegetation coverage, thus has low soil and water conservation function. On the other hand, the physical and chemical properties of soil are easily affected by human activities, so the combination of human activities and soil factors has a greater impact on vegetation.

Considering the results of factor detector and interactive detector, we conclude that the physical factors, rather than the anthropogenic factors have been the primary driving factors for vegetation dynamics across northern Tibet during 2000 and 2015. This result consists with the previous findings conducted in the similar ecosystems [10,42]. The anthropogenic factors, especially population density and grazing density, also play important roles in vegetation variation. In our study, precipitation, temperature, and sunshine duration are the dominant physical factors influencing vegetation variation. Especially, the decreasing trend of precipitation, increasing trend of temperature, and decreasing trend that sunshine duration may cause an unsuitable environment for vegetation growth. In this case, we indicate that the effects of the anthropogenic factors may be magnified when climate conditions become unfavorable for vegetation. As we all know, the ecosystem in the Qinghai-Tibet plateau are vulnerable and extremely sensitive to the changing environment. We can not regulate climate in the short term, but corresponding measures and adaptive strategies should be taken to regulate the negative impacts of human activities, and ultimately to achieve sustainable development of the alpine ecosystems.

5. Conclusions

In summary, based on the NDVI time series, this study explored the issue of vegetation variation on the alpine grassland ecosystem. Our study clearly showed a significant decreasing trend of vegetation index across northern Tibet during the period of 2000–2015. The physical factors, rather than the anthropogenic factors, have been the primary driving factors for vegetation dynamics, while the effects of the anthropogenic factors may be magnified when the physical environment is unsuitable for vegetation growth. Specifically, meteorological factors were the primary factors that determined vegetation variation in the alpine meadow and the alpine scrub, and the anthropogenic factors influenced the alpine scrub variation the most. The interactions between the driving factors had a larger effect on vegetation than any single factor. In the alpine meadow, the areas with high altitude may improve the impact of precipitation on vegetation. Meanwhile, decreasing precipitation combined with increasing temperature and decreasing sunshine duration may potentially enhance the degradation of the grasslands, while the soil types and the anthropogenic factors together have resulted in a reduction of grassland plant cover in the alpine steppe. In addition, the influence of population density is highly similar to that of grazing density, which may provide access to simplify the study of the anthropogenic activities in the plateau ecosystem. Although this study integrally quantified the influence of environment factors on vegetation variation, it failed to identity the dominant factors at a sub-region scale. Finer resolution data are needed for further studies.