Spatiotemporal Differentiation and the Factors of Ecological Vulnerability in the Toutun River Basin Based on Remote Sensing Data

Ecological vulnerability assessment increases the knowledge of ecological status and contributes to formulating local plans of sustainable development. A methodology based on remote sensing data and spatial principal component analysis was introduced to discuss ecological vulnerability in the Toutun River Basin (TRB). Exploratory spatial data analysis and a geo-detector were employed to evaluate the spatial and temporal distribution characteristics of ecological vulnerability and detect the driving factors. Four results were presented: (1) During 2003 and 2017, the average values of humidity, greenness, and heat in TRB increased by 49.71%, 11.63%, and 6.51% respectively, and the average values of dryness decreased by 165.24%. However, the extreme differences in greenness, dryness, and heat tended to be obvious. (2) The study area was mainly dominated by a high and extreme vulnerability grade, and the ecological vulnerability grades showed the distribution pattern that the northern desert area was more vulnerable than the central artificial oasis, and the central artificial oasis was more vulnerable than the southern mountainous area. (3) Ecological vulnerability in TRB showed significant spatial autocorrelation characteristics, and the trend was enhanced. The spatial distribution of hot/cold spots presented the characteristics of “hot spot—cold spot—secondary hot spot—cold spot” from north to south. (4) The explanatory power of each factor of ecological vulnerability was temperature (0.5955) > land use (0.5701) > precipitation (0.5289) > elevation (0.4879) > slope (0.3660) > administrative division (0.1541). The interactions of any two factors showed a non-linear strengthening effect, among which, land use type ∩ elevation (0.7899), land use type ∩ precipitation (0.7867), and land use type ∩ temperature (0.7791) were the significant interaction for ecological vulnerability. Overall, remote sensing data contribute to realizing a quick and objective evaluation of ecological vulnerability and provide valuable information for decision making concerning ecology management and region development.


Introduction
In recent years, the ecosystem is undergoing increasing pressures and degeneration owing to climate change, urban sprawl, and human activities [1][2][3][4]. The environmental problems following those pressures, such as global warming [5], soil erosion [6], desertification [7], environment pollution [8], and loss of biodiversity [9], have changed the ecosystem structure and process, decreased the ecosystem service value, and posed great threats to region sustainable development. Ecological vulnerability is regarded as a core issue in sustainability science [10], and it is an inherent property of an ecosystem. As a Sustainability 2019, 11, 4160 3 of 19 differences moisture index (NDMI) from Landsat 5 TM data to identify the water content pixel and its vulnerability level in the Keleghai River Basin, India. Hu et al. [38] established an integrated index to evaluate ecological vulnerability in Fuzhou City based on data from Landsat ETM+/OLI/TIRS while the impacts of human activities on ecosystems inevitably spread to deserts and mountains. However, the application of remote sensing data in ecological vulnerability evaluation is still at an exploration stage, and the situation needs to be improved. Most studies have demonstrated vulnerability with the land use type and landscape index from remote sensing images but have paid less attention to vulnerability as an inherent attribute or state of ecosystem from a raster scale. An arid region is a vast and sparsely populated area where precipitation is scarce, the temperature is high, and humidity is low [39]. These climate characteristics and other biophysical features have profound implications on the ecological environment and society economy and also easily result in critical ecological environmental decline and long-term vulnerability that is harmful to ecosystem stability and human life [40][41][42]. Especially in the arid land of northwestern China, human activities are concentrated in an artificial oasis, from where, from just a small percentage of arid oasis, while the impacts of human activities on ecosystems inevitably spread to deserts and mountains. Until now, few studies have attempted to adequately assess ecological vulnerability and spatial heterogeneity from the perspective of assessing ecological status in a raster scale.
The Toutun River originates from the Kalawucheng Mountain and flows about 190 km northward into Junggar Basin. The river flows through Urumqi City, Wujiaqu City, and Changji City, containing a typical mountain-oasis-desert ecosystem in an arid area. The Toutun River Basin (TRB) is a typical arid oasis in Xinjiang, northwest China, with a fragile ecological environment, while it possesses an important position in regional development and is strategically significant due to its position on the country's periphery [43]. In 2017, the population and gross national product accounted for 11.83% and 28.51% of Xinjiang, respectively, while the land area accounted for 1.39% of Xinjiang. However, because of the typical temperate continental arid climate and increasing development intensity in the middle oasis, urban expansion [43], lack of water resources [44], environment pollution [45], soil erosion [46], desertification [47], and the fragile ecological environment have already negatively influenced the local development and posed great challenges to the sustainable development of the basin, making the TRB appropriate for assessing ecological vulnerability. Establishing a quantitative ecological vulnerability assessment framework in a raster scale is important for a more detailed understanding of the ecosystem in TRB. Research on spatial heterogeneity and its driving mechanism will help in gaining a deep and systematic understanding of the distribution, changes, and interactions in ecological vulnerability. These will enrich the evaluation methods of ecological vulnerability and promote the remote sensing data application in ecological vulnerability evaluation. In this study, we used the humidity index, greenness index, dryness index, and heat index extracted from remote sensing data to evaluate ecosystem vulnerability with the spatial principal component analysis method. Then, exploratory spatial data analysis and a geo-detector model were employed for further study of ecological vulnerability. Specific objectives were included to (1) assess ecological vulnerability with remote sensing data in TRB; (2) analyze the temporal and spatial distribution characteristics of ecological vulnerability; (3) judge the explaining power of influence factors and the interaction between factors.

Study Area
TRB is located in the center of the Xinjiang Uygur Autonomous Region (86 • 24 32" N-88 • 09 53" N, 43 • 56" E-45 • 20 01" E), with an area of 23,090.62 km 2 (Figure 1). Backed by the Tianshan Mountain, it lies on the southern edge of the Gurbantunggut Desert with an elevation from 274 to 4512 m. The elevation gradually increases from north to south. The area consists of a desert and an oasis and is mountainous, which is typical in arid inland river basins. In the temperate continental arid climate zone, the study area has a short spring and autumn, a cold winter, a hot and dry summer, and little precipitation all year round, what is typical in arid/semiarid areas in northwest China. The area has strong solar radiation, less precipitation and high evaporation, and a distinct soil-vegetation vertical belt formed under the influence of the climate and topography.

Data Source and Pre-Processing
A database was established including remote sensing data, digital elevation model (DEM) data, meteorological data, and land use data for the study of the ecological vulnerability in TRB. Thematic layers ware established with the same spatial boundary, a WGS_1984_UTM_45N projected coordinate system, and a resolution of 30 m. (1) Remote sensing data. Landsat images with a resolution of 30 m in 2003, 2010, and 2017 were obtained from the United States Geological Survey. After that, the remote sensing images were preprocessed by ENVI 5.3 software produced by Esri. Using a radiometric correction tool, the digital number of the original image was switched to reflectance at a sensor, and FLAASH atmospheric correction reduced the deviations caused by light and the atmosphere. (2) DEM data and meteorological data. Both of them were downloaded from the Resources and Environment Data Cloud Platform (http://www.resdc.cn/). Then, the elevation and slope were extracted from DEM data. (3) Land use data. In accordance with the land use classification system in China and the characteristics in TRB, land use type was divided into the following categories: Desert, water bodies,

Data Source and Pre-Processing
A database was established including remote sensing data, digital elevation model (DEM) data, meteorological data, and land use data for the study of the ecological vulnerability in TRB. Thematic layers ware established with the same spatial boundary, a WGS_1984_UTM_45N projected coordinate system, and a resolution of 30 m. (1) Remote sensing data. Landsat images with a resolution of 30 m in 2003, 2010, and 2017 were obtained from the United States Geological Survey. After that, the remote sensing images were preprocessed by ENVI 5.3 software produced by Esri. Using a radiometric correction tool, the digital number of the original image was switched to reflectance at a sensor, and FLAASH atmospheric correction reduced the deviations caused by light and the atmosphere. (2) DEM data and meteorological data. Both of them were downloaded from the Resources and Environment Data Cloud Platform (http://www.resdc.cn/). Then, the elevation and slope were extracted from DEM data. (3) Land use data. In accordance with the land use classification system in China and the characteristics in TRB, land use type was divided into the following categories: Desert, water bodies, forest, grassland, cropland, construction land, bare exposed rock or gravel, and glacier and perennial snowfield. Regarding the GoogleEarth image and field survey data, the overall precisions of the land use/cover results were above 85%, which satisfied the research needs.

Ecological Vulnerability Indicators Retrieval
In arid areas, water condition and vegetation cover are the main factors affecting the characteristics of the ecological environment. Moisture, greenness, dryness, and heat are widely applied to assess ecological status due to their close relationship with ecology quality. The ecological vulnerability index (EVI) is constructed to assess ecological vulnerability with remote sensing data based on remote sensing ecological index (RSEI) proposed by Xu [48].

• Wetness retrieval
The Tasseled Cap transformation is widespread when monitoring and assessing ecosystems [49]. The components of Tasseled Cap, including brightness, greenness, and wetness, have an immediate impact on the physical parameters of land surface [50]. The wetness component is exploited for land surface moisture (LSM) here, and the formulas based on the reflectance of Landsat TM and Landsat OLI image are expressed as follows, respectively: where b Blue , b Green , b Red , b NIR , b SWIRI , and b SWIR2 is the planetary reflectance of the blue, green, red, near-infrared, and mid-infrared band of the Landsat 5 TM image or Landsat 8 OLI image, respectively.

• Greeness retrieval
Vegetation growth is a description of status that reflects the confrontation between natural ecosystems and drought [24]. In general, the near-infrared and red bands of remote sensing images are recommended to quantitatively reflect the growth of photosynthetic vegetation on the ground [51,52]. As an effective indicator for quantifying vegetation, the Normalized Difference Vegetation Index (NDVI) is capable of reflecting greenness and calculated as follows: where b Red and b NIR denote the planetary reflectance of the red and near-infrared band of the Landsat 5 TM image and Landsat 8 OLI image, respectively.

• Dryness retrieval
The Impervious Built-up Index (IBI) and Soil Index (SI) can effectively distinguish built-up lands and bare areas, respectively. In consideration of urban construction land, deserts, and abandoned land in the study area, a combination index of IBI and SI, named Normalized Difference Imperviousness and Soil Index (NDISI), is applied due to it being adept at accurately recognizing impervious ground and bare ground without masking out water body [49,53].
is the planetary reflectance of the blue, green, red, near-infrared, and mid-infrared band of the Landsat 5 TM image and Landsat 8 OLI image, respectively.
• Heat retrieval The thermal infrared band of remote sensing data contains information on land surface temperature [54]. A significant correlation between satellite-derived surface temperature and the actual land surface temperature has been proven in many studies of urban heat island [55,56]. The thermal Band 6 of Landsat 5 TM and the thermal Band 10 of Landsat 8 OLI were employed to calculate the land surface temperature (LST) in 2003, 2010, and 2017. The calculations of LST are expressed as follows [57]: L λ = gain * DN + bias where L λ is the at-satellite spectral radiance values of the thermal band; gain is the gain value of the thermal infrared band, and bias is the offset value, both of which can be read from the header file of the corresponding image. T b is the at-satellite brightness temperature, and K 1 and K 2 are the thermal conversion constants; at the thermal infrared band of Landsat 5, K 1 = 607.76 W/(m 2 .sr.µm) and K 2 = 1260.56 K; at the thermal infrared band of Landsat 8, K 1 = 774.89 W/(m 2 .sr.µm) and K 2 = 1321.08 K. λ is the wavelength of the thermal infrared band; ρ = 1.4380 * 10 4 µm; ε is the surface specific emissivity.

Construction of EVI
To remove the influence of different dimensions and ranges of values from the initial variables, original data were normalized in a range from 0 to 1 for "homogenization" and "nondimensionalization". If there was a positive correlation between indicators and EVI, the equation was expressed as follows: While, if there was a negative correlation between indicators and EVI, the equation was: where X i is the normalized value of variable i, and x i is the initial value of variable i.
Principal component analysis can reduce the dataset dimensionality by converting the raw data attributes into linearly uncorrelated variables [58]. Statistical methods have always brought much inconvenience to the study of spatial data. Considering the obstacle in processing spatial and temporal variations, a combination of the geographic information system and principal component analysis, defined as spatial principal component analysis [59], which is adept at processing spatial data, was employed here. The principal components were chosen based on the following criterion: If the accumulated contribution of eigenvalues controlled 85% of the total, the principal components were reserved, and the rest of the components could be ignored. A composite index EVI for an objective and quantitative evaluation is expressed as follows: where r i is the contribution ratio of principal component (PC), and i is the quantity of PC that remains.
The contribution ratio r i is calculated as follows: where p i represents the eigenvalue of PC i . The evaluation models of the ecological vulnerability index in 2003, 2010, 2017 were obtained from Table 1: The distribution of ecological vulnerability from 2003 to 2017 was obtained by the raster calculation tool in Arcgis10.5. Combined with the characteristics of the TRB, ecological vulnerability values were classified into five categories with natural breaks as follows: Grade I-no vulnerability, Grade II-slight vulnerability, Grade III-moderate vulnerability, Grade IV-high vulnerability, and Grade V-extreme vulnerability. The Ecological Vulnerability Body Index (EVBI) was used for analyzing the overall differences in ecological vulnerability, the formula is: where P i is the ecological vulnerability grade; no vulnerability, slight vulnerability, moderate vulnerability, high vulnerability, and extreme vulnerability are assigned to 1, 2, 3, 4, and 5, respectively; A i represents the area of grade i; S represents the total area of TRB.

Exploratory Spatial Data Analysis
As an integration of methods for describing and visualizing the spatial distribution of observations, exploratory spatial data analysis describes spatial distributions and heterogeneity by identifying the spatial outliers or discovering spatial correlation patterns [60]. In practical application, Global Spatial Autocorrelation and Hot Spot (Getis-Ord Gi*) Analysis are often applied to explore the spatial distribution characteristics of observation data. • Global Spatial Autocorrelation Global Spatial Autocorrelation measures if and how much the dataset is autocorrelated throughout the study region. Global Moran's I is calculated as follows [61,62]: x i ; x i and x j represent the values of spatial units i and j, respectively; n is the number of spatial units; w ij is an element of the spatial weight matrix W, if spatial units i and j share a common border, w ij = 1; otherwise, w ij = 0. The value of Global Moran's I is in the range of −1 to 1.
• Hot spot (Getis-Ord Gi*) Analysis Getis-Ord Gi* is a useful technique for describing the spatial cluster of high values and low values [60,63]. The G * i statistic is calculated as follows: x j , n, and w ij are the same as above. The G * i values were divided into four grades with natural breaks from high to low. After that, four categories were proposed as follows: Hot spot area, secondary hot spot area, secondary cold spot area, and cold spot area. If high values cluster, G * i will be large, and it will be defined as a hot spot, while if low values cluster, G * i will be small, and it will be defined as a cold spot.

Geo-Detector Model
Geo-detectors are a set of statistical methods for detecting spatial differentiation among the geographical elements [64]. This paper employed factor detection and interaction detection to analyze the explanatory power and the interaction between multiple variables. The explanatory power of variable X on attribute Y was calculated with the following formula: where q is the explaining power of a variable on a spatial attribute; t is the categories or partitions of variables; N h is the quantity of sample units in subfields; n is the quantity of sample units in the whole area; L is the quantity of subfields; σ 2 is the variance of a single variable of the entire region; and σ 2 h is the variance of the subfield.
Interaction judgement makes a contribution by distinguishing between the interactions between any two variables. Interaction detection was employed here to detect whether the various factors formed a combined influence on ecological vulnerability, or whether the effect between different variables was independent. The basis for interaction type is shown in Table 2.

Characteristics of Indicators
The statistics of four indicators are shown in Table 4.  From the spatial distribution of ecological factors (Figure 2), it can be seen that in the desert areas in the north, the soil moisture content decreased, and the surface temperature increased. Because of the expansion of cultivated land at the edge of the middle artificial oasis, the ecological factors were all developing in a favorable direction. The vegetation cover and soil moisture in the south mountainous area were improved, and the tendency to expand to high altitude areas was obvious. Sustainability 2019, 11, x FOR PEER REVIEW 10 of 18

Ecological Vulnerability Status
The ecological vulnerability of the basin has been reduced as a whole; the EVBI in 2003, 2010, and 2017 were 3.6337, 3.5031, and 3.5393, respectively. According to Figure 3, ecological vulnerability grades showed the distribution pattern that the northern desert area was more vulnerable than the central artificial oasis, and the central artificial oasis was more vulnerable than the southern mountainous area. In 2003-2017, the vulnerability level was dominated by high vulnerability and extreme vulnerability, mainly distributed in the northern desert area and gradually decreasing with

Ecological Vulnerability Status
The ecological vulnerability of the basin has been reduced as a whole; the EVBI in 2003, 2010, and 2017 were 3.6337, 3.5031, and 3.5393, respectively. According to Figure 3, ecological vulnerability grades showed the distribution pattern that the northern desert area was more vulnerable than the central artificial oasis, and the central artificial oasis was more vulnerable than the southern mountainous area.
In 2003-2017, the vulnerability level was dominated by high vulnerability and extreme vulnerability, mainly distributed in the northern desert area and gradually decreasing with the expansion of the artificial oasis. The area surrounding the artificial oasis was dominated by moderate vulnerability and slight vulnerability. The artificial oasis development increased the efficiency of water; the desert landscape was gradually replaced by the cropland landscape, and then, the area with no vulnerability and slight vulnerability gradually expanded along the edge of the artificial oasis. On the other hand, the encroachment of construction land on cropland strengthened the vulnerability in the interior of the oasis. The ecological vulnerability grades of the southern mountainous areas were mainly no vulnerability and slight vulnerability. As a result of global warming, the rising temperature accelerated ablation of glaciers, increased the river runoff, and vegetation was limited. The ecological vulnerability in southern mountainous areas was also weakened. the expansion of the artificial oasis. The area surrounding the artificial oasis was dominated by moderate vulnerability and slight vulnerability. The artificial oasis development increased the efficiency of water; the desert landscape was gradually replaced by the cropland landscape, and then, the area with no vulnerability and slight vulnerability gradually expanded along the edge of the artificial oasis. On the other hand, the encroachment of construction land on cropland strengthened the vulnerability in the interior of the oasis. The ecological vulnerability grades of the southern mountainous areas were mainly no vulnerability and slight vulnerability. As a result of global warming, the rising temperature accelerated ablation of glaciers, increased the river runoff, and vegetation was limited. The ecological vulnerability in southern mountainous areas was also weakened. The area of each ecological vulnerability grade was calculated and is shown in Figure 4. In 2003, the area with no vulnerability (I), slight vulnerability (II), moderate vulnerability (III), high vulnerability (IV), and extreme vulnerability (V) was 1566.77 km 2 , 2546.89 km 2 , 3070.64 km 2 , 3321.98 km 2 , and 6599.42 km 2 , respectively. From 2003 to 2010, the area with no vulnerability, slight vulnerability, and high vulnerability increased to 1887.41 km 2 , 2713.60 km 2, and 3949.20 km 2 , respectively, while the area with moderate vulnerability and extreme vulnerability increased to 2982.91 km 2 and 5572.57 km 2 , respectively. Area conversion mainly occurred between adjacent levels; for example, the increasing area of no vulnerability, slight vulnerability, and high vulnerability came from slight vulnerability, moderate vulnerability and extreme vulnerability, and high vulnerability, respectively. From 2010 to 2017, the area with no vulnerability, slight vulnerability, and extreme vulnerability increased to 2002.23 km 2 , 3060.90 km 2 and 6547.97 km 2 , respectively, while the area with moderate vulnerability and high vulnerability decreased to 2300.62 km 2 and 3193.97 km 2 , respectively. The main types of conversion included high vulnerability to extreme vulnerability, moderate vulnerability to slight vulnerability, and slight vulnerability to no vulnerability. The area of each ecological vulnerability grade was calculated and is shown in Figure 4. In 2003, the area with no vulnerability (I), slight vulnerability (II), moderate vulnerability (III), high vulnerability (IV), and extreme vulnerability (V) was 1566.77 km 2 , 2546.89 km 2 , 3070.64 km 2 , 3321.98 km 2 , and 6599.42 km 2 , respectively. From 2003 to 2010, the area with no vulnerability, slight vulnerability, and high vulnerability increased to 1887.41 km 2 , 2713.60 km 2, and 3949.20 km 2 , respectively, while the area with moderate vulnerability and extreme vulnerability increased to 2982.91 km 2 and 5572.57 km 2 , respectively. Area conversion mainly occurred between adjacent levels; for example, the increasing area of no vulnerability, slight vulnerability, and high vulnerability came from slight vulnerability, moderate vulnerability and extreme vulnerability, and high vulnerability, respectively. From 2010 to 2017, the area with no vulnerability, slight vulnerability, and extreme vulnerability increased to 2002.23 km 2 , 3060.90 km 2 and 6547.97 km 2 , respectively, while the area with moderate vulnerability and high vulnerability decreased to 2300.62 km 2 and 3193.97 km 2 , respectively. The main types of conversion included high vulnerability to extreme vulnerability, moderate vulnerability to slight vulnerability, and slight vulnerability to no vulnerability.
respectively. From 2010 to 2017, the area with no vulnerability, slight vulnerability, and extreme vulnerability increased to 2002.23 km 2 , 3060.90 km 2 and 6547.97 km 2 , respectively, while the area with moderate vulnerability and high vulnerability decreased to 2300.62 km 2 and 3193.97 km 2 , respectively. The main types of conversion included high vulnerability to extreme vulnerability, moderate vulnerability to slight vulnerability, and slight vulnerability to no vulnerability.

Spatial Heterogeneity in Ecological Vulnerability
Ecological vulnerability index maps were used to analyze the Global Moran's I of ecological vulnerability. As shown in Figure 5, the Global Moran's I in 2003, 2010, and 2017 passed the significance test, implying that ecological vulnerability in TRB had a positive autocorrelation or highly clustered pattern. From the perspective of the evolution trend, the spatial autocorrelation degree was continuously strengthened, and the spatial autocorrelation level was enhanced.

Spatial Heterogeneity in Ecological Vulnerability
Ecological vulnerability index maps were used to analyze the Global Moran's I of ecological vulnerability. As shown in Figure 5, the Global Moran's I in 2003, 2010, and 2017 passed the significance test, implying that ecological vulnerability in TRB had a positive autocorrelation or highly clustered pattern. From the perspective of the evolution trend, the spatial autocorrelation degree was continuously strengthened, and the spatial autocorrelation level was enhanced. G * was classified into four categories with natural breaks, and four categories were proposed as follows: Hot spot areas, secondary hot spot areas, secondary cold spot areas, and cold spot areas ( Figure 6). The distribution of cold/hot spots in TRB was relatively stable and presented strong spatial accumulation. On the whole, the cold/hot spots of the ecological vulnerability of the TRB generally showed the characteristics of "hot spot-cold spot-secondary hot spot-cold spot" from north to south. The distribution of cold spots was mainly concentrated in the south of TRB, while hot spots were distributed in the north; meanwhile, significant differences were presented in spatial distribution patterns and quantities of cold and hot spots. G * i was classified into four categories with natural breaks, and four categories were proposed as follows: Hot spot areas, secondary hot spot areas, secondary cold spot areas, and cold spot areas ( Figure 6). The distribution of cold/hot spots in TRB was relatively stable and presented strong spatial accumulation. On the whole, the cold/hot spots of the ecological vulnerability of the TRB generally showed the characteristics of "hot spot-cold spot-secondary hot spot-cold spot" from north to south. The distribution of cold spots was mainly concentrated in the south of TRB, while hot spots were distributed in the north; meanwhile, significant differences were presented in spatial distribution patterns and quantities of cold and hot spots.
The distribution of cold/hot spots showed typical zonal features. The hots pots and secondary hot areas were mainly distributed in the desert and desert steppe between the oasis and mountain, and their spatial structure presented a blocky distribution. The cold spots and secondary cold points were widely distributed in the mountainous areas in the south and east of TRB, where there is always low ecological vulnerability owing to the high altitude, abundant precipitation, and flourish vegetation. In the artificial oasis, agricultural land and urban green land would be beneficial to the ecological environment, while building land, forming a dry and exposed surface, would have an adverse impact on the ecological environment. Because of this and the above-mentioned factors, cold/hot spots of vulnerability in the artificial oasis were distributed in a complex, point-like distribution structure.  G * was classified into four categories with natural breaks, and four categories were proposed as follows: Hot spot areas, secondary hot spot areas, secondary cold spot areas, and cold spot areas ( Figure 6). The distribution of cold/hot spots in TRB was relatively stable and presented strong spatial accumulation. On the whole, the cold/hot spots of the ecological vulnerability of the TRB generally showed the characteristics of "hot spot-cold spot-secondary hot spot-cold spot" from north to south. The distribution of cold spots was mainly concentrated in the south of TRB, while hot spots were distributed in the north; meanwhile, significant differences were presented in spatial distribution patterns and quantities of cold and hot spots. The distribution of cold/hot spots showed typical zonal features. The hots pots and secondary hot areas were mainly distributed in the desert and desert steppe between the oasis and mountain, and their spatial structure presented a blocky distribution. The cold spots and secondary cold points were widely distributed in the mountainous areas in the south and east of TRB, where there is always low ecological vulnerability owing to the high altitude, abundant precipitation, and flourish

Relationship Between Detected Factors and Ecological Vulnerability
Based on the results of the geographical detector (Table 5), the q values of factors were presented as follows: Temperature (0.5955) > land use (0.5701) > precipitation (0.5289) > elevation (0.4879) > slope (0.3660) > administrative division (0.1541), and all factors were statistically significant with respect to the changes in the ecological vulnerability. As a typical basin oasis in the arid area, the TRB was sensitive to hydro-thermal combination. Appropriate temperature and abundant rainfall played an important role in ecological vulnerability by increasing vegetation cover and inhibiting the surface temperature. The land use type reflected human activities to a certain extent. Grassland, woodland, and cultivated land were covered by vegetation, the humidity and greenness of which were higher than those of desert and build-up land. The "oasis effect" and "cold island" effects also had a significant effect on the oasis ecosystem's response to arid climates. Due to the combined use of water and thermal resources, artificial oases had a high production efficiency and showed a stable ecological state. Concerning the significance test, the temperature, land use, precipitation, elevation, slope, and administration division were selected to detect the interaction effect in ecological vulnerability. The interactions between any two indicators of ecological vulnerability were higher than the explanatory power of each single indicator but lower than the sum of two indicators' explanatory power ( Table 6). The result indicated that the interaction of any two factors in ecological vulnerability showed non-linear strengthening. The maximum value of interaction appeared between land use type and elevation, which meant that the difference of ecological vulnerability was the largest in different elevation under the same land use type or different land use types at the same elevation. The land use combined with other factors significantly enhanced the explanatory power of the factors. Among them, the interaction of land use type ∩ elevation (0.7899), land use type ∩ precipitation (0.7867), and land use type ∩ temperature (0.7791) were the significant control factors for ecological vulnerability in the TRB.

Discussion
Redundancy of multispectral data is often encountered and needs to be properly handled [65]. The combined index EVI was selected to reflect dryness, heat, wetness, and greenness. At the same time, the correlation and redundancy between data was reduced by spatial principal component analysis. The presence of multicollinearity can be diagnosed by examining its indicators. The variance inflation factor (VIF) and tolerance index (TOL), defined as complementary to VIF, were employed to estimate and remove potentially redundant features. The cutoff value of VIF is 10 [66]; if VIF > 10 (i.e., TOL < 0.1), the multicollinearity among indicators are severe [66]. Taking the data of 2017 as an example, the VIF and TOL of indicators were calculated with 17,089 points related to EVI and indicators, and then, the results showed that the VIF was less than 10, and the TOL was more than 0.1, which indicated that there was no obvious correlation between the indicators, and the indicators and methods were desirable (Table 7). Ecological vulnerability assessment can objectively reveal the spatial differentiation features of ecological vulnerability in TRB and provide a scientific basis and references for the ecological environmental protection management of the basin. This method presented an apparent advantage compared with other studies and can be applied to assess vulnerability, particularly in arid lands. In TRB, ecological vulnerability showed a distribution pattern which not only corresponded with the actual situation, a typical mountain-oasis-desert pattern in an arid area, but was also consistent with the development requirements of the region [44]. Adequate and rational desert construction in arid areas expanded the living space of human beings and improved the quality of the local ecological environment [67,68]. With the development of water and soil resources at the edge of the artificial oasis, the desert landscape in alluvial plain was replaced with an artificial oasis, and the ecological environment improved along the edge of the artificial oasis. Reservoir construction and inter-basin water transfer projects have effectively solved the problem of uneven water resources distribution; however, the redistribution of water also means changes to the original ecological pattern. In arid areas where water resources are relatively scarce, over-exploitation of land (mainly including over-expansion of cultivated land and cities) means more water is being consumed for agricultural production and city development. As a consequence, the natural oasis ecosystems may not receive the needed water, and desertification, salinization, groundwater level reduction, and vegetation cover reduction occur easily [69][70][71]. Additionally, with the cities' extension, croplands became occupied by construction land, while the negative impacts of urban sprawl such as a decrease in biodiversity, greenhouse gas emissions, and air pollution also followed. From the long-term view, over-exploitation of cultivated land and cities will threaten oasis stability and lead to oasis development being unsustainable. Given the complexity of the ecosystem, some microscopic indicators and ecological processes, such as socio-economic impacts, vegetation community structure, and soil degradation also need to be focused on. In future research on ecological vulnerability, we would consider more attributes of the ecosystem in arid areas and verify the result with statistical data and field survey data.

Conclusions
Ecological vulnerability assessment is of great significance for ensuring and improving ecosystem stability. This paper made an attempt to evaluate ecological vulnerability using only remote sensing data at a raster scale, which should be useful for a rapid and objective understanding of ecological status. In this paper, NDVI, LSM, LST, and NDISI were retrieved from remote sensing images, and then, EVI was constructed with spatial principal component analysis to obtain the distribution characteristics of the ecological vulnerability in TRB. The results showed that the average values of humidity, greenness, and heat increased, and the average values of dryness decreased. However, the extreme differences in greenness, dryness, and temperature tended to be obvious from 2003 to 2017. Although the ecological vulnerability of TRB has been reduced as a whole, it was mainly dominated by high and extreme grades of ecological vulnerability. The ecological vulnerability grades showed the distribution pattern that the northern desert area was more vulnerable than the central artificial oasis, and the central artificial oasis was more vulnerable than the southern mountainous area. Ecological vulnerability had significant spatial autocorrelation characteristics, and the trend was enhanced. The distribution of cold/hot spots showed typical zonal features and presented the characteristics of "hot spot-cold spot-secondary hot spot-cold spot" from north to south. The explanation for each factor of ecological vulnerability was temperature > land use > precipitation > elevation > slope > administrative division. The interactions of any two factors presented a non-linear strengthening effect. Among them, the interaction of land use type ∩ elevation, land use type ∩ precipitation, and land use type ∩ temperature were the significant control factors for ecological vulnerability.
The framework presented in this paper provides a visual and mensurable approach for a detailed understanding of ecological vulnerability in a raster scale, and the method can be employed to other areas with suitable datasets for ecological conservation and environmental management, particularly in arid regions. Furthermore, field survey data and statistical data will provide a reference for ecological vulnerability assessment and verify the accuracy of assessment; they should be taken into account in future research.