A Long-Term Ecological Vulnerability Analysis of the Tibetan Region of Natural Conditions and Ecological Protection Programs

: The combined impacts of drastic natural environment change and increasing human interference are making the uncertainty of the Tibetan Plateau’s ecological vulnerability the world’s largest. In this study, an ecological vulnerability index ( EVI ) of Tibet in the core area of the Tibetan Plateau was assessed using a selected set of ecological, social, and economic indicators and using a spatial principal component analysis (SPCA) to calculate their weights. The data included Landsat images and socio-economic data from 1990 to 2015 in ﬁve-year intervals. The results showed that the total EVI remained at a high vulnerability level, with drastic ﬂuctuation from 1990 to 2000 (a peak in 1995, when there was a sudden increase in light vulnerability, which moved to extreme vulnerability in the next period), and minor ﬂuctuations after 2000, gradually increasing from southeast to northwest. In addition, the spatial analysis showed a distinct positive correlation between the EVI and grassland area (0.33), land use degree (0.15), NDVI (0.14), livestock husbandry output, and a negative correlation in terms of desertiﬁcation area. The artiﬁcial afforestation program (AAP) had a positive signiﬁcant correlation with NDVI (R 2 = 0.88), preventing the environment from becoming more vulnerable. The results provide practical information and suggestions for planners to improve the land use degree in urban areas and the vegetation coverage in pastoral regions of the Tibetan Plateau based on the spatial–temporal heterogeneity patterns of the EVI of Tibet.


Introduction
The world has been focusing on keeping the ecological ecosystem in a relatively stable dynamic equilibrium state in different aspects, such as the need to balance economic development and ecological conservation [1], and the ecosystem and communities' response to exterior pressure [2,3]. When the external influences irreversibly exceed the balance range, environmental problems occur worldwide, such as severe climate anomalies [4,5], land desertification [6], ecosystem degradation [7], and natural disasters [8], especially in ecologically vulnerable regions. The potential of an ecosystem's capability to respond, cope with, and recover when the ecosystem is disturbed by a specific hazard or stressor in a specific time and space is called ecological vulnerability [9][10][11][12][13][14]. Vulnerability includes exposure, sensitivity, and adaptive capacity, which vary according to human interaction with nature. An ecological vulnerability assessment is one of the ways to assess climate change, a highlight in ecological research in recent years [15][16][17]. Therefore, by identifying and targeting vulnerable regions, the threat of intense external forces can be mitigated or prevented, demonstrating the crucial importance of assessing ecological vulnerability for sustainable management and development.
The Intergovernmental Panel on Climate Change (IPCC) presented an effective assessment structure to evaluate ecological vulnerability accurately in 2001 [18]. Ecological vulnerability is considered as the inherent property of the ecosystem to react to hazard sensitivity and its resilience to cope with unexpected pressure [19][20][21]. Furthermore, the methods for assessing ecological vulnerability have been modified based on this structure in the 21st century. Typically, evaluating ecological vulnerability is based on three parts: sensitivity, resilience, and pressure (the PSR model) [22][23][24]. Studies in Europe and coastal areas have mainly focused on the sensitivity of the ecosystem and its response to climate change and pattern change over time [25][26][27][28][29][30], which measure ecosystem tolerance and vegetation response predictions to climate anomalies. In mainland areas, research on ecological vulnerability is more concentrated on land use patterns [31][32][33], desertification soil problems [24,[34][35][36], and human disturbances [37][38][39][40][41]. Typically, vulnerable zones with different topological characteristics are studied [38][39][40] to determine the determinants of vulnerability. Appropriate ecological management policy should consider both environmental and socio-economic factors, which are built on an objective and integrated evaluation system. Therefore, a comprehensive ecological vulnerability evaluation system with high accuracy should be applied to these vulnerable areas [41][42][43][44][45]. Scientists have made progress in evaluation methods, including scenario analysis methods [46], quantitative evaluation model methods, and the most widely used, comprehensive index methods. The comprehensive index evaluation system contains selected indicators and the weights of each indicator. The studies predominantly use the analytic hierarchy process (AHP) [24,47,48] and expert scoring [48] to calculate the weight. However, these methods are not just non-objective but also controversial and may not apply to other areas. The analytical unit matters significantly as it displays the accuracy of the study. Yet, currently, most studies use an administration unit instead of a unitive spatial grid unit, resulting in rough conclusions with high inaccuracy and a lack of spatial comparison [49][50][51][52].
As the most fragile and vital region in the world, the Tibetan Plateau, also known as the third pole, has experienced the severe threat of a damaged ecological security barrier and has witnessed serious ecological problems including freeze-thaw erosion [53], hydraulic erosion [54], land desertification [34], grassland degradation [3,55], and salinization [56]. These issues have made the ecosystem difficult to inhabit and maintain a healthy state. The Tibet Autonomous Region is at the margin of the Tibetan Plateau. Approximately 92% of the entire autonomous region lies within the arctic-alpine environment [3] and suffers from frequent crustal movements, drastic wind erosion [8], and freezing erosion [53]. The melting of snow and ice is the fundamental prerequisite for freeze-thaw erosion, and solar radiation is the decisive factor for melting snow and the doubling of the global average in warming speed [51]. Grasslands cover more than 65% of the total Tibet land area [3], and, due to the alpine drought environment, the grassland ecosystem has undergone a significant degradation, leading to decreases in vegetation diversity and productivity [55,57]. Since the grassland ecosystem is irreplaceable for ecological security and pasture output, the inference from climate change and human activities (forage grazing and urbanizing) becomes exceptionally crucial. Hence, the government has implemented several ecological protection projects to increase vegetation coverage in Tibet (Table 1), and the artificial afforestation program (AAP) is one of the most important and effective means of land greening, especially in arid areas [34,47], which has increased the vegetable coverage rate of Tibet by 16.9% [58].
The AAP is largely implemented in all ecological protection projects listed in Table 1 including country-wide programs such as the Grain for Green Program, the National Comprehensive Demonstration Area for Sand Control and Prevention, and the Closing Hill for Afforestation Program, etc. These listed programs were launched by the Chinese government in the 1980s and have been continuously invested in for the past 25 years [34,47,58]. For instance, the Protection and Construction of Ecological Security Barrier Plan in Tibet, which was conducted by the Chinese government in 2008, had a total investment from the government of 56. 66  are still disputes over the use of AAP, as the cost efficiency of AAP is not ideal, and the operation period of AAP exceeds the expectation, which has dismissed the actual needs of local pastoralists [59,60]. Therefore, more investigation on the specific effects and extent of AAP's influence on reducing ecological vulnerability in Tibet is needed. Mainly, current ecological vulnerability studies in Tibet are limited to specific prefectures or certain aspects of driving factors such as the Sichuan-Tibet railway [37]. These studies have inaccurate analytical units, a small range of study periods, influences from single-type indicators, and inadequate evidence on determinant analysis [40][41][42][43][44][45]. Though numerous studies have investigated ecological health and ecosecurity and assessed the ecological risk of the Tibetan Plateau, more research is still needed on an ecological vulnerability analysis with an ecological protection project policy in the Tibetan Plateau [1,15,31,35,45,49].
In this study, the method was integrated to establish a multiobjective multidimensional ecological vulnerability evaluation system and to analyze the driving factors and relevance between ecological vulnerability and ecological protection project policies. This paper aims to (1) assess EVI in spatial-temporal scale and its determinants in Tibet accurately using spatial principal components analysis (SPCA), (2) identify the relationship between EVI and grasslands and urbanization, (3) validate the impacts that AAP brings to EVI and provide valuable and practical information for planners in environmental policymaking in the Tibetan Plateau.

Study Area
Tibet autonomous region (Figure 1), located in western China and the west and south of the Qinghai-Tibetan Plateau (26 • 50 N-36 • 53 N, 78 • 25 E-99 • 06 E), accounting for more than half of the area of the Tibetan Plateau and about one-eighth of the total area of China (1.23 × 10 6 km 2 ), known as the "roof of the world" and the "third pole of the earth", is the highest altitude place in the world [3]. The region's topography can be divided into three major areas: the northern Tibetan plateau, the Brahmaputra River basin, and the eastern Tibetan canyon area [57,61]. These regions provide Tibet with abundant natural resources. In these fertile places, the light resources, wind resources, and water resources are the most significant, which are twice than the amount of east coast China in 2016, the average number of windy days per year is up to 100-150 days, with a maximum of up to 200 days, and the largest theoretical reserves of hydro energy resources in China amount to 2.1 × 10 8 KW, respectively [62,63].

Evaluation System Framework
The comprehensive ecological vulnerability evaluation system combines the study area's ecological and socioeconomic conditions in a modified pressure-state-response structure [19][20][21][43][44][45]57,64]. Furthermore, the evaluation index system was constructed by referring to the index system of ecological vulnerability evaluation in the Tibetan high land mountains pre-established by related scholars and expert consultants familiar with the local ecological environment [1,6,15,21,57,62]. The evaluation system is established on the basis of the accessibility and the continuity of the data [64][65][66][67]. The methodologica flowchart is shown in Figure 2.

Evaluation System Framework
The comprehensive ecological vulnerability evaluation system combines the study area's ecological and socioeconomic conditions in a modified pressure-state-response structure [19][20][21][43][44][45]57,64]. Furthermore, the evaluation index system was constructed by referring to the index system of ecological vulnerability evaluation in the Tibetan highland mountains pre-established by related scholars and expert consultants familiar with the local ecological environment [1,6,15,21,57,62]. The evaluation system is established on the basis of the accessibility and the continuity of the data [64][65][66][67]. The methodological flowchart is shown in Figure 2.

Evaluation System Framework
The comprehensive ecological vulnerability evaluation system combines th area's ecological and socioeconomic conditions in a modified pressure-state-re structure [19][20][21][43][44][45]57,64]. Furthermore, the evaluation index system was cons by referring to the index system of ecological vulnerability evaluation in the Tibeta land mountains pre-established by related scholars and expert consultants famili the local ecological environment [1,6,15,21,57,62]. The evaluation system is establis the basis of the accessibility and the continuity of the data [64][65][66][67]. The methodo flowchart is shown in Figure 2.
As vulnerability is largely interpreted through mountain hazards, grassland degradation, and human activities, variables were chosen with a predisposition to these hazardaffected elements. For example, there are water erosion problems on the southeast edge of Tibet [54] where the environmental factors related to water erosion are topography, vegetation, and precipitation, and the corresponding indicators are slope, vegetation type, vegetation cover, and precipitation. The indices reflected the extent to which these environmental problems affected local ecological vulnerability. The environmental factors related to wind erosion are topography, vegetation condition, and wind speed, and the corresponding indicators are roughness, vegetation cover, and average wind speed.
The evaluation system included 18 variables. These indicators were selected from four major impact factors, which were climate, topography, vegetation, and human disturbance (slope, topographic relief, vegetation coverage, land use degree, landscape diversity, desertification area, plateau permafrost area, grassland area, water resource area, average annual precipitation, relative humidity, average annual temperature, hours of sunshine, wind speed, solar radiation intensity, population density, GDP, and livestock husbandry output). There was information redundancy in these different aspects' interactions, which was eliminated using a spatial principal components analysis.
This comprehensive ecological vulnerability evaluation system demonstrates the influence of socioeconomic factors on ecological vulnerability, the extent of each indicator, and its relationship with ecological vulnerability.
Considering the influence of different grading schemes on the final vulnerability evaluation results, and the strong subjectivity of grading schemes, to avoid excessive interference of human factors in the process of index grading, this paper adopts the evaluation model of normalized weighted summation of fuzzy classification of indicators for comprehensive research. Continuous data, such as vegetation cover (0, 1), are normalized by the extreme value method, and all indicators with different scales and ranges of values are unified to (0, 1). The effect directions express whether each indicator has a positive or negative relationship with ecological vulnerability. The positive indicator uses the value of the normalization, and the negative indicator is the value of 1 minus the normalization value.
The negative indicators include NDVI, landscape diversity, grassland area, water resource area, average annual precipitation, relative humidity, and average annual temperature. The positive indicators include slope, topographic relief, land use degree, desertification area, plateau permafrost area, hours of sunshine, solar radiation intensity, population density, GDP, and livestock husbandry output. Land use degree represents the proportion of developed land area including built-up area, farmland area, and forest area among the total area in one spatial grid.
Therefore, all indicators are standardized to (0, 1) before the overlay analysis with a 2 km spatial grid evaluation unit [64].
Moreover, in Tibet, the current ecological vulnerability evaluation research is primarily based on the county evaluation unit instead of a uniform spatial grid unit [42], resulting in rough conclusions with high inaccuracy that lack the necessary detailed spatial analysis and localization.

Weight Calculation
Climate, geography, vegetation, social development, and economics are ecological impact elements that may interact in this study's multi-index integrated ecological vulnerability assessment methodology. Due to the exchanges, there is information redundancy. In the analysis process, this would result in overlapping and redundant information, as well as erroneous and problematic analytical findings [10,11,15,16,21,68,69]. The SPCA technique may successfully maintain relevant primary components by removing information overlap in the original indicator data [12,21,36,43]. The indicators' weights are determined using SPCA (Spatial Principal Component Analysis) based on PCA (principal component analysis) after the standardization of all indicators in ArcGIS 10.4 [64][65][66][67][68]. Its calculation principle is consistent with PCA, which transforms multiple variables into a few principal components through dimensionality reduction and produces a correlation coefficient matrix composed of each standardized index [10][11][12]44,60,[66][67][68][69][70][71][72].
A total of 18 standardized indices were analyzed by principal components analysis to generate 4 new comprehensive indices. In addition, the absolute value of eigenvalue presents the extent of the correlation of indicators and principal components; the more significant the eigenvalue is, the stronger the relationship will be [11,12,17,36,68,72]. The first four principal components whose cumulative contribution rate reached 85% or more were selected; the final principal component result is shown in Table 2, and the whole contribution of original indexes of ecological vulnerability to principal components is shown in Appendix A, Table A1. The calculation formula was as follows: where R was the correlation coefficient matrix, Z was the standardized value of each selected index, n was the number of indexes, λ was the eigenvalues of the R correlation coefficient matrix, I was the identity matrix, CCR was the cumulative contribution rate, m was the number of principal components that were determined, P was the matrix containing values of every considered principal component, and W was m number of eigenvectors with the largest eigenvalues selected to form the matrix.

EVI Calculation and Gradation
According to the above method, the higher the EVI, the more vulnerable the ecological environment. The index is obtained by the sum of the comprehensive values of multiple principal components and their corresponding weights [11,12,36,[60][61][62][63][64][68][69][70], as shown in the following formula: where EVI is the ecological vulnerability index; r is the contribution ratio; P is the principal component; m is the number of principal components; r i is the contribution ratio of the i principal component; and n i is the eigenvalue of the i principal component.
To quantitatively assess the changing trend of the vulnerability of regional eco-systems, a comprehensive regional ecological vulnerability index (EVSI, Ecological Vulnerability Standard Index) needs to be constructed and standardized according to Equation (7) to obtain the comprehensive ecological vulnerability index EVSI of Tibet:

Hot Spot Analysis of Ecological Vulnerability
Hot Spot analysis is a spatial analysis model used to display the degree of spatial clustering calculated by Getis-Ord Gi*. The results of the EVI calculations are visualized using the cold/hot spot model. In this study, the Getis-Ord Gi* index was used to analyze the high/low spatial aggregation degree of EVI changes, that is, the spatial distribution of cold/hot spots. Moreover, cold/hot spot analysis is the database that underpins heterogeneity analysis. The calculation formula was: where Gi* was the output statistical Z-score, x j was the EVI change of space unit j, and w ij was the spatial weight between adjacent space units i and j.

Spatial Correlation Analysis between EVI and Urbanization Index
According to a previous study [13,41], the composite urbanization index (UI) can be quantified through population urbanization, economic urbanization, and land urbanization. These three levels are expressed via population density (PD), GDP density (GDPD), and the area occupied by built-up land (ULP), respectively.
where UI represents the urbanization index of the evaluation unit and PD , GDPD , and ULP represent the population density, GDP density, and construction land area ratio of the evaluation unit after standardization. The spatial correlation between EVSI and UI was analyzed using Moran's I index. The Monte Carlo test was applied to the analysis 999 times. If Moran's I > 0, it means that there is a positive correlation between them, and vice versa, there is a negative correlation; additionally, the spatial clustering types of EVSI and UI were obtained in GeoDa software, which were classified as insignificant: high-high, low-low, low-high, and high-low. The formula is as follows: where I is the global bivariate Moran's I for the EVSI and urbanization and I i kl is the local bivariate Moran's I for the EVSI and urbanization level. N stands for the total number of spatial units. W ij stands for spatial weight matrix for measuring the spatial correlation between the i and j spatial units. Z i refers to the deviation between the attribute of i spatial unit and the average of the attribute. Z j refers to the deviation between the attribute of j spatial unit and the average of the attribute. X i k refers to the value of attribute k of spatial unit i. X k refers to the average of attribute k. σ k is the variance of attribute k. X j l refers to the value of attribute l of spatial unit j. X l refers to the average of attribute l. σ l refers to the variance of attribute l.

Correlation between NDVI and Afforestation Area and Regression Tests
The NDVI reflects the absorption and wavelengths characteristics of vegetation in the red and near-infrared regions and, therefore, provides a good indication of ground vegetation growth [73]. The annual and seasonal NDVI averages reflect vegetation growth during the year and season. NDVI is an essential indicator in Tibet, especially when evaluating the grassland growing status and forest conservation [73][74][75]. In addition, the grassland area is the most significant driving factor in ecological vulnerability, according to the result in Section 3.3. Statistical analysis methods must be applied to figure out the effectiveness of AAP policy and the relationship between the effect and EVI. Therefore, NDVI can also reflect the effectiveness of the AAP policy and the correlation between the AAP policy and EVI changing patterns. The NDVI can be the indicator to reflect the extent of the ecosystem coping with the exposure to a hazard [76,77]. However, it is not only determined by the internal resilience of the system, but also by the anthropogenic improvement such as afforestation programs [34,47]. During the study period, some early afforestation area increases was not consistent with NDVI, while the increases in afforestation area alongside decreases in NDVI (especially in 2000, 2005). The regression models in Statistical Product and Service Solutions (SPSS) software were used in this study to analyze the correlation type to explore the relationship between NDVI and afforestation area using the continuous data from 1990 to 2017. The time series is considered a linear regression function of time: Y = at + b (t stands for the year; where a is the tendency value, a > 0 indicates an increasing trend with time; a < 0 indicates a decreasing trend with time, and the absolute value of a value reflects the rate of increase or decrease.).

Spatial and Temporal Changes of Ecological Vulnerability
The EVI was classified into five levels using the natural breakpoint classification method (Jenks) in ArcGIS (slight vulnerability: <1.0527, light vulnerability: 1.0527-1.3672, medium vulnerability: 1.3672-1.6817, high vulnerability: 1.6817-2.0676, and extreme vulnerability: >2.0676). Throughout the study period, the EVI of Tibet values varied from 0.423 to 2.418, with its highest value in 1995 at 1.925, lowest in 1990 at 1.573, and mean value at 1.774, indicating that Tibet was in a high vulnerability state as a whole.
The proportion of areas with different types of ecological vulnerability for 1990-2015 is shown in Figure 3. On the temporal scale, the light vulnerability area proportion had a short peak in 1995 at 47%, with the most prominent area having an increasing fluctuation of 258.53%. However, the area's proportion of light vulnerability soon decreased to 9% in 2000. The extreme vulnerability area presented an increasing pattern during the study period and accounted for the most significant area proportion from 2000 onward at around 44%. The area distribution ratio of five ecological vulnerability levels stayed stable in descending order of extreme vulnerability, medium vulnerability, high vulnerability, slight vulnerability, and light vulnerability.   Table 3 shows the percentages of fluctuation of areas with different EVI valu administrative prefecture. The Ngari region was the most vulnerable in all six p and Nagqu experienced a peak in 2005 at 2.088, equal to the EVI of Ngari. In 2000 than half of the prefectures presented an EVI decrease, including Lhasa, Qamdo, Sh and Nyingchi (the most significant one). The EVI in 1995 was significantly reduce  Table 3 shows the percentages of fluctuation of areas with different EVI values per administrative prefecture. The Ngari region was the most vulnerable in all six periods, and Nagqu experienced a peak in 2005 at 2.088, equal to the EVI of Ngari. In 2000, more than half of the prefectures presented an EVI decrease, including Lhasa, Qamdo, Shannan, and Nyingchi (the most significant one). The EVI in 1995 was significantly reduced, and all went through a decline. The spatial distribution of continuous EVI is shown in Figure 4. Generally, the EVI of Tibet increased primarily in the northern area (mainly in Nagqu and Ngari) over 25 years. Furthermore, the EVI presented a progressively downward trend from north to south. Ecological vulnerability deteriorated considerably in Ngari, Nagqu, and Xigaze. Additionally, the EVI of Ngari stayed in an extremely high status during the study period (high vulnerability in 1990, extreme vulnerability in 1995-2015). Shannan, Lhasa, and Xigaze experienced a significant increase over 25 years, and Xigaze changed from light vulnerability to medium vulnerability. The area of the Brahmaputra River basin had an evident increment in EVI with a spatial pattern of downward trending along the west-east axis, indicating the lower reaches of Brahmaputra River (mainly Nyingchi City) were steadily low vulnerable. The light and slight vulnerability areas were mainly distributed in valley networks and the adret slope of the Transhimalaya and Himalaya, namely, the relatively low-altitude areas.

Spatial Heterogeneity Analysis of Ecological Vulnerability
The distribution of cold/hot spots of ecological vulnerability in Tibet during the study period is shown in Figure 5.
Overall, the hot spot was in a stable status with a north to south descending trend. The cold spot generally decreased, mainly in the southern Xigaze, northern Shannan, and Lhasa. The EVI of the whole Nyingchi was in a 99% confidence cold state, indicating that the null hypothesis of complete spatial randomness should be rejected, and this area was 99% expected to present a significant low EVI clustered pattern. Temporal changes demonstrated that the hot spot showed a decreasing trend from 1990 to 2015, with the largest area of hot spots reaching its peak in 1990. The cold spots showed a decreasingincreasing trend from 1990 to 2015, with the largest area of cold spots in 1990. Generally, the hot spot occupied a larger area in Tibet during the study period. Meanwhile, the area of the hot spot and cold spot showed a slight fluctuation after 1995, indicating increasing ecological capacity and stability in cold spot areas.
Regions with significant EVI change filled around 73% of the total Tibet area. Moreover, these regions were mainly located on the edge of Tibet on the north-south axis. Northern Xigaze, Southern Nagqu, and Northern Qamdo showed stable status over 25 years with no significant EVI changes. Figure 5 shows that the hot spot areas were mainly concentrated in Ngari and Nagqu, and the cold spot areas in Nyingchi, Shannan, and Southern Qamdo. During the study periods, the cold spot in Xigaze, Lhasa, and Shannan was decreased, with an overall EVI increment in these prefectures. In the basin of Brahmaputra, the cold spot area diminished significantly from the middle part to the east part, while the hot spot area stayed unfluctuating.

Spatial Heterogeneity Analysis of Ecological Vulnerability
The distribution of cold/hot spots of ecological vulnerability in Tibet during the study period is shown in Figure 5.

Determinant Indicators of EVI
The principal component analysis results are shown in Figure 6. The first four principal component layers, which contributed over 88%, were used to explore the drivers of change in ecological vulnerability in Tibet. Table 1 lists the primary environmental protection programs constructed in Tibet during the study years, which may be relevant to the indices influencing the EVI in Tibet.
In principal component 1 (PC1), the EVI was significantly positively correlated with land use factors such as grassland area (0.3484), NDVI (0.2978), and meteorology factors such as average annual precipitation (0.2983), and average annual temperature (0.2983) from 1990 to 2015. The economic factors are relatively negative with the EVI, with livestock output (−0.1440). Moreover, solar radiation intensity (−0.1379) and wind speed (−0.2135) played the most negative role in 25 years, especially the wind speed factor in 1990 of −0.2530. Since the PC1 accounted for more than half the proportion of the final principal component (55.522%), the indicators that affected the PC1 also played a significant role in the integrated driving factors. Grassland area, average annual precipitation, and average annual temperature were the most vital factors correlated with EVI in Tibet.
PC2, PC3, and PC4 of the EVI were positively correlated with livestock output, land use degree, desertification area, wind speed, and solar radiation intensity, which were quite different from PC1. In general, the grassland area factor remained the most evident one to influence the EVI. The land use degree factor became critical after 1990; the desertification area played an important role, especially in 2010 and 2015. Specifically, GDP was negatively connected to EVI in 1990, while the correlation became diametrically opposite in from 1995 to 2015.
Analysis of the integrated correlation index of PC1-PC4 in 1990-2015 showed that the most influential factor of the EVI was the grassland area, especially in 2005 and 2010. This was followed by land use degree, livestock output, desertification area, and NDVI. Overall, the hot spot was in a stable status with a north to south descending trend. The cold spot generally decreased, mainly in the southern Xigaze, northern Shannan, and Lhasa. The EVI of the whole Nyingchi was in a 99% confidence cold state, indicating that the null hypothesis of complete spatial randomness should be rejected, and this area was 99% expected to present a significant low EVI clustered pattern. Temporal changes demonstrated that the hot spot showed a decreasing trend from 1990 to 2015, with the largest area of hot spots reaching its peak in 1990. The cold spots showed a decreasingincreasing trend from 1990 to 2015, with the largest area of cold spots in 1990. Generally, the hot spot occupied a larger area in Tibet during the study period. Meanwhile, the area

Effects of Urbanization on EVI
From 1990 to 2015, both the EVI and the urbanization level in Tibet increased significantly, indicating a synergistic change over timescales. The results showed that the EVI and the urbanization level were positively correlated in 1990-2015, and the positive correlation gradually strengthened. The Z-values in Table 4 showed that the EVI and urbanization level had a powerful spatial aggregation effect, with the weakest only in 1990.

Determinant Indicators of EVI
The principal component analysis results are shown in Figure 6. The first four principal component layers, which contributed over 88%, were used to explore the drivers of change in ecological vulnerability in Tibet. Table 1 lists the primary environmental protection programs constructed in Tibet during the study years, which may be relevant to the indices influencing the EVI in Tibet.    Table 5 shows the urbanization index value of all prefectures in Tibet from 1990 to 2015. The indices barely fluctuated in 1995-2015 with a nearly ±0.001 swing. Qamdo, Nagqu, and Ngari were the only prefectures that had grown in urbanization level since 1990. In contrast, Nyingchi had the most significant drop in 1995 from 0.481 (1990) to 0.296 (1995). Table 5 indicates that the urbanization level had spatial heterogeneity, and even though the urbanization level was positively connected to the EVI, the situation was different in each prefecture.

Correlation between NDVI and Afforestation Area
The interannual trend line equations of the NDVI and afforestation area are listed below in Table 6, indicating that NDVI and afforestation area were positive from 1990 to 2015. Fisher's test was applied, and the result showed a significant relationship between the NDVI and afforestation area with the R 2 > 0.85, p < 0.01, and F > 19.00. The standard values of EVI, NDVI, and afforestation area from 1990 to 2015 are shown in Figure 7.

Correlation between NDVI and Afforestation Area
The interannual trend line equations of the NDVI and afforestation area are lis below in Table 6, indicating that NDVI and afforestation area were positive from 1990 2015. Fisher's test was applied, and the result showed a significant relationship betw the NDVI and afforestation area with the R 2 > 0.85, p < 0.01, and F > 19.00. The stand values of EVI, NDVI, and afforestation area from 1990 to 2015 are shown in Figure 7.   Figure 7 show that the NDVI in Tibet over the 25 years had a posit trend in all, and the afforestation area had an obviously positive trend. This trend me a better turn in the vegetation growth status in afforested areas, and the spatial distri tion density also showed an increasing trend in the last 25 years. The negative releva   Figure 7 show that the NDVI in Tibet over the 25 years had a positive trend in all, and the afforestation area had an obviously positive trend. This trend means a better turn in the vegetation growth status in afforested areas, and the spatial distribution density also showed an increasing trend in the last 25 years. The negative relevance between the afforestation area and the EVI is shown in Figure 7, indicating that the EVI decreased every five years when the afforestation area increased.

The Spatial-Temporal Patterns of Ecological Vulnerability
Overall, the EVI in Tibet has been in a high vulnerability state for 25 years, while extreme vulnerability since 2000 occupied more than 43% of the total study area. The results identified the dynamic patterns in time series and spatial heterogeneity of EVI in Tibet.
In terms of time series, the proportion structure of different ecological vulnerability levels was unstable until 2000. Many high vulnerability areas turned into light vulnerability, and part of them changed into extreme vulnerability. This sudden change in 1995 showed the highest ecological vulnerability value in Tibet at 1.925 throughout the study period. The hot spot analysis showed that in 1995 the count of the not significant area increased and reached a peak (81,153 pixels) with nearly the same amount of cold spot (86,839 pixels), indicating that the increment of ecological vulnerability was mostly spatially random except the north corner of Tibet. Accordingly, the Ngari prefecture and Nagqu prefecture, located in the north part of Tibet, experienced a sharp increase in EVI in 1995. In the meantime, the slightest ecological vulnerability lay in southeastern Tibet, with Nyingchi having the smallest EVI. With Ngari possessing the highest EVI in Tibet from 1990 to 2015, Ngari was considered the most vulnerable region in Tibet, which is in accordance with the hyperaridity, low precipitation, slight hypsography physiognomy topography fluctuation, and severe glacier erosion ecosystems. Moreover, the Ali Shiquanhe Township Phase I and Phase II Sand Control and Sand Management Project was conducted from 1992 to 2000 in Ngari, resulting in urbanization development and relevant construction [78]. At the same time, the construction made the urbanization index increase accordingly (Table 7) and decrease light vulnerability areas, which made the environment more vulnerable by disarraying resource development. In the meantime, according to Table 3, southeastern Tibet had a better ecosystem status in comparison to northwestern Tibet, with Shannan and Nyingchi having the least EVI in Tibet during 1990-2015. Though the EVI in southern Tibet rose suddenly in 1995 with the increasing rate of Shannan (58%), Lhasa (45%), Xigaze (36%), it was still less than the northern part and generally stayed in a light vulnerability status. In southeastern Tibet, though the threat of aridity and the significant difference in temperature and elevation, permanent permafrost was diminished, there were also other issues that deteriorated the environment. For instance, abundant forest resources led to over woodcutting and forest damage. Since the timber for the whole of Tibet is mostly taken from southeastern Tibet, this doubled the fall amount in comparison to the growth of forest trees. The situation of unbalanced tree felling brought an increase in the desertification rate and therefore increased the EVI of southern Tibet. There was an evident aggregation in river valleys, low-elevation areas, the adret slope of mountains, and water bodies regarding the slight and light vulnerability area. This clustering certified that protecting these areas is essential. Moreover, the government approved a series of policies to enhance the ecological security barrier of Tibet since 2000, which made the ecosystems in Tibet more stable, yet still fragile. These policies mainly concentrated on solving the desertification issue and improving vegetation cover rate through building conservation areas and expanding the grassland and forest area.
The downstream of the Brahmaputra River basin and open-water bodies were stable and had a low vulnerability, while the upper reaches were the opposite. This consequence means that the Brahmaputra River basin demonstrated a significant spatial variation pattern with EVI descending by the river flow. These findings are consistent with previous research.
Overall, the variation of EVI polarized from 1990 to 2015, with northern Tibet having more fragile in zones with low landscape diversity and arid alpine situations.

Driving Factors of Ecological Vulnerability
The ecological vulnerability driving factors can be traced back to natural factors and human interferences. As the result of correlation analysis between principal components and various indicators shown in Figure 6, grassland area, land use degree, desertification area, livestock output, and NDVI are the top five critical indicators that affect ecological vulnerability. Solar radiation intensity and wind speed account for the least weight. This proportion indicated that the green vegetation prosperity was the dominant fragment in all factors, followed by the human disturbance indicators (land use degree and livestock output). In comparison, the significance of climate factors and topological factors were relatively inessential.
Among all selected indicators, the influence of the grassland area factor was large not only because the grassland area accounted for the most significant proportion in Tibet, but because of the ecosystem functions of grassland in conserving water and soil, maintaining biodiversity, and controlling wind and sand preventing the environment from becoming more vulnerable [55,57]. The grassland in northwestern Tibet is different from the southeastern. Alpine meadow wetlands form most ecosystems in northwestern Tibet, surrounded by bare lands and glaciers, water bodies, and forests barely distributed there. Therefore, this idiosyncratic landscape pattern may lead to a dramatic response to external interferences and irreversible vulnerability deterioration.
Moreover, as the environment is austere in northwestern Tibet, its population density is relatively low in accordance with an artificial disturbance that placed more responsibilities on the natural aspect resulting in the highest ecological vulnerability index. Meteorological factors including precipitation, relative humidity, temperature, wind speed, and hours of sunshine were the fundamental ingredients affecting the growth of vegetation dynamics [79]. This fact is then related to soil erosion caused by rainfall, wind, and melted snow [80]. Combined soil erosion with topological elements such as elevation and slope can cause brutal geological disasters [81], which further aggravates the vulnerability of the area. However, the situation is different in southeastern Tibet. Even though the desertification problem still exists here, the cause was not similar to north Tibet. Abundant heat and rainfall along with low elevations and water resources bring south Tibet various landscape patterns and high vegetation coverage. The ecosystem here is far more habitable than in the north, with the essential standards for living guaranteed. Consequently, the urbanization level was lifted by pastoral area intensification and massive project constructions such as the Sichuan-Tibet Railway, diminishing the area of slight and light vulnerability. The disorderly series of constructions raised the urbanization rate of Tibet, according to Table 5. Relatively, EVI was lifted along with the urbanization rate; as a result, Section 3.4 showed a continuous significant positive relation between urbanization rate and EVI. Remarkably, this pattern did not fit every prefecture. In Shannan, the correlation between urbanization and EVI was negative during 1990-2015 [82], which meant that though EVI in Tibet overall was influenced by urbanization, the situation was spatially heterogeneous.
Furthermore, the correlation was gradually strengthened, meaning the influence of urbanization on EVI was fortified and needed more attention to balancing urbanized construction and ecosystem protection. In addition to the disturbance of urbanization, desertification influence was inescapable. The desertification area gradually increased from 1990 to 2015 and was prominently distributed in river valleys and adret regions with lower elevation and closer distance to rivers. This situation appeared in the lower Brahmaputra River basin and valleys resulting from water swilling, geological disasters, and artificial project constructions, especially hydropower projects.
The government has been taking a series of systematic measures to diminish human activities' influence and prevent natural factors from deteriorating the environment (Table 1). These policies mainly concentrated on increasing vegetation coverage by building or recovering forest and grassland, especially artificial afforestation construction. The positive correlation between the Artificial Afforestation Programme and NDVI to a further extent indicated that AAP raised the vegetation coverage and, therefore, ameliorated the ecosystem status. The improving trend demonstrated critical evidence that with the conduction of AAP, EVI changes and forest areas became better, and this connection strengthened over 25 years, emphasizing the effectiveness of the AAP.
In general, EVI in Tibet is largely influenced by grassland area, and also by land use degree and NDVI. Moreover, land use degree had more significant influence in southeastern Tibet than northwestern Tibet, as the urbanization rate in southeastern Tibet exceeded that of northwestern Tibet.

Sustainable Implications for Ecosystems Management
The objective reality of Tibet's fragile ecological environment and the pressure of economic development have made environmental protection in Tibet a severe challenge. As the core variable of environmental change in the Pan-Third Pole, the Tibetan Plateau will pose a significant challenge to the survival and development of more than 20 countries and three billion people in the Pan-Third Pole region if adequate measures are not taken [60]. This study showed that the EVI of Tibet is significantly negatively affected by grassland areas, which are influenced by extreme climate and massive project construction. Therefore, this article offers the following suggestions to prevent the ecological vulnerability of Tibet from increasing further.
Firstly, the grassland ecosystem is susceptible to climate change and grazing. Soil erosion in northern Tibet should be paid attention to in order to take precautions against mountain hazards. These places have been in a highly vulnerable state for 25 years, and even though the government has approved ecological policies to take protection governance, the result still showed a worsening trend. This trend means more powerful and practical measures are needed. Grazing bans and a grass-livestock balance system should be conducted strictly and improving the pasture contract management system can boost the artificial grass planting and natural grassland improvement project and development of shed-feeding and semi-shed-feeding farming. In terms of extreme vulnerability areas in northern and northwestern Tibet, grassland ecological protection subsidy and incentive mechanisms and returning pasture to grass would reinforce the current ecological protection. For southern Tibet, the seriously degraded grasslands need to be fenced off from grazing permanently, and moderately degraded grasslands need to be fenced off from grazing temporarily. Moreover, as the results of artificial afforestation have shown a positive effect on the ecosystem, this policy should be continued and strengthened to accelerate the promotion of sandy land management in middle Tibet and northwest Tibet, and vigorously prevent and control soil erosion in important slight and light vulnerability areas (i.e., water-conserving areas and water and wind erosion staggered areas).
Secondly, from the perspective of optimizing land use degree, pastoral management, refined land function planning, and change of industry structure are needed to solve the aggravated fragmentation of the ecological system. The central conflict here is the need for exploitation and restrictive ecosystem protection. Therefore, in those fenced areas, the primary function can be to plant crops to ensure meeting food demand while the environment is guarded. Moreover, expanding green areas, waters, and wetlands can build a balanced habitable living system as the ecological space is increased. These operational measures can optimize the land use degree and, therefore, improve the ecological vulnerability accordingly.
In summary, the government has implemented a series of environmental protection policies over the past 25 years, and most of them have achieved good results, especially the desertification control along the Brahmaputra River. The government should maintain the above ecological and environmental protection achievements while exploring more effective and appropriate methods for local development.

Conclusions
Ecological vulnerability in Tibet varies greatly both temporally and spatially with changes in the natural landscape fragmentation and anthropogenic disturbances such as grazing. This study analyzed the temporal and spatial variation, heterogeneity, driving factors, and the positive aspects that environmental policies have contributed to EVI through a new, integrated evaluation system to quantify ecological vulnerability. The analysis leads to the following conclusions: (1) Tibet's ecological vulnerability level experienced a drastic fluctuation from 1990 to 2000 (peaking in 1995); the proportion pattern stabilized since 2000 with a descending value order as extreme vulnerability (occupied nearly half of Tibet), medium vulnerability, high vulnerability, slight vulnerability, and light vulnerability. In general, Tibet was in a high vulnerability state in the study period. The EVI is more spatially than temporally embodied. The spatial distribution of EVI mainly manifested as a higher vulnerability state in the northwest than that in the southeast with a decreasing trend.
In the Brahmaputra River Basin, the EVI increased over 25 years with a descending order, upstream to downstream. (2) During the study period, the hot spot center stayed in north Tibet, but the significance of clustering increased; in the meantime, the cold spots clustered in the southeast part of the study area in 1990-2015. The hot spots and the cold spots reached the largest area in 1990, indicating gradual stabilization of the EVI after 1995. This study provides a scientific basis for natural and land use resource development, ecological environment protection, and planning and construction. Therefore, this study can provide policy implementation recommendations for coordinating ecosystem protection and socioeconomic development, helping Tibet to achieve a sustainable state and improve its ecological vulnerability to benefit the Qinghai-Tibet Plateau.