**2. Materials and Methods**

### *2.1. Overview of the Study Area*

The Hexi Corridor, which has a length of 900 km and width of 100 km (Figure 1), is located in the arid region of northwestern China. It is bordered by Wushaoling to the east, the Jade Gate to the west, the Qilian and Altun mountains to the south, and the Mazongshan, Heli, and Longshou mountains to the north. Owing to its temperate continental climate, annual precipitation is in the range of 40–300 mm, and annual temperature is in the range of 6.2–9.0 ◦C. It should be noted that annual evaporation exceeds 1500 mm in most areas. The land use types in this region include farmland, forest, grassland, desert, water bodies, snow or glaciers, and residential land.

### *2.2. Data Sources*

This study used 30-m spatial resolution Landsat5/TM data acquired in 2000, Landsat8/OLI data acquired in 2020, and MOD13A3 data acquired during 2000–2020, which were sourced from both the Aerospace Information Research Institute of the Chinese Academy of Sciences (www.aircas.cas.cn, accessed on 1 February 2021) and the National Aeronautics and Space Administration (www.nasa.gov, accessed on 1 February 2021). Additionally, 1-km spatial resolution PANDA data acquired during 2000–2020 were provided by the National Tibetan Plateau Data Center (http://data.tpdc.ac.cn, accessed on 1 May 2021).

**Figure 1.** Location of meteorological stations and DEM information in the study area.

To analyze the regional climate change, daily temperature and precipitation data recorded during 2000–2020, at 19 meteorological stations in the Hexi Corridor, were obtained from the Gansu Meteorological Bureau (http://cdc.cma.gov.cn/, accessed on 1 February 2021).

Digital elevation model data with a 30-m spatial resolution (ASTER GDEM V3) were sourced from the Geospatial Data Cloud (http://www.gscloud.cn/, accessed on 1 May 2021), which is a product of jointly developed by Japan's METI and US NASA.

### *2.3. Research Methods*

The data processing and analysis flow are displayed in Figure 2.

### 2.3.1. Analysis of Land Use Types

Classification of the land use in the Hexi Corridor was performed on the basis of the Landsat5/TM, Landsat8/OLI, and digital elevation model data. The original images were first preprocessed, by means of orthophoto correction, geometric correction, radiometric calibration, atmospheric correction, mosaic and clipping, cloud removal, and shadow processing. Generally, machine-learning approaches, hybrid methods, fuzzy theory, and other methods were the means of obtaining land use information. On the basis of the landscape classification system, supervised classification and manual correction were used to categorize the land use in the Hexi Corridor during 2000–2020 into the following types: farmland, forest, grassland, residential land, bare land or desert, water bodies, and snow or glaciers. Three hundred random samples were extracted to validate the classification of two different images, whose interpretation matched reality, with total accuracies and Kappa coefficients of 87.67% and 0.857 in 2000, and 89.33% and 0.872 in 2020.

A transition matrix can depict in detail the spatial structure and transformation direction of different land use types for different time sequences. This facilitates an accurate understanding of the direction of change of the type, structure, and distribution ratio of land use [15]. The land use classification of the Hexi Corridor during 2000–2020 underwent data fusion, and then an intersection analysis was performed to obtain the transition data for each land use type, to create a land type transition matrix of the overlapping areas.

**Figure 2.** Data processing and analysis flow chart.

### 2.3.2. Index of Landscape Pattern Indexes

To investigate the changes of landscape pattern over the past 21 years, number of patches (NP), mean patch area (AREA), landscape shape index (LSI), aggregation index (AI), Shannon's diversity index (SHDI), and contagion index (CONTAG) were calculated [16].

### 2.3.3. NDVI Change Trend Analysis

The maximum value composite method was used to maximize the bimonthly synthetic standard vegetation index, and to obtain annual normalized difference vegetation indexes (NDVIs). Then, following removal of zero and negative values, the annual NDVIs for different vegetation types were obtained by overlaying vegetation maps.

A linear regression analysis was conducted, to simulate the spatial variation of the annual NDVIs for the different vegetation types in the target area [17,18]. The calculation formula can be expressed as follows:

$$\theta\_{slope} = \frac{n \times \sum\_{j=1}^{n} j \times NDVI\_j - \sum\_{j=1}^{n} j \sum\_{j=1}^{n} NDVI}{n \times \sum\_{j=1}^{n} j^2 \left(\sum\_{j=1}^{n} j\right)^2} \tag{1}$$

where *n* refers to the cumulative number of monitored years, *NDVIj* is the NDVI per year in year *j*, and *θslope* signifies the slope of the trend line. A value of *θslope* >0(*θslope* < 0) indicates that the trend of the NDVI change over *n* years is increasing (decreasing).

The Theil–Sen Median and Mann–Kendall test were combined to perform a long-term sequence analysis of vegetation change [19–21]. The Mann–Kendall test has been effectively applied in temporal sequence analyses in the fields of hydrology and meteorology. The Theil–Sen median can be used to assess whether temporal sequence data show an upward or downward trend [22,23]. The median of the slope in *n* (*n* − 1)/2 data combinations was calculated in the trend analysis. The calculation formula can be expressed as follows:

$$S\_{NDVI} = Median\,\,\frac{NDVI\_j - NDVI\_i}{j - i},\tag{2}$$

where 2000 ≤ *i* < *j* ≤ 2020. A value of *SNDVI* >0(*SNDVI* < 0) indicates that the NDVI is increasing (decreasing).

When an NDVI is tested using the Mann–Kendall test, the NDVI of a certain temporal sequence is regarded as a set of independently distributed sample data, and parameter *Z* serves as the attenuation index of the NDVI pixel. The calculation formula can be expressed as follows:

$$Z = \begin{cases} \frac{S-1}{\sqrt{s(S)}}, & S > 0 \\ & 0, \ S = 0 \\ \frac{S+1}{\sqrt{s(S)}}, & S < 0 \end{cases} \tag{3}$$

$$S = \sum\_{j=1}^{n^\*-1} \sum\_{i=j+1}^n \text{Sgn}\left(\text{NDVI}\_{\bar{j}} - \text{NDVI}\_i\right),\tag{4}$$

$$\text{Sign}(\text{NDVI}\_{\text{j}} - \text{NDVI}\_{\text{i}}) = \begin{cases} 1, & \text{NDVI}\_{\text{j}} - \text{NDVI}\_{\text{i}} > 0 \\ 0, & \text{NDVI}\_{\text{j}} - \text{NDVI}\_{\text{i}} = 0, \\ -1, & \text{NDVI}\_{\text{j}} - \text{NDVI}\_{\text{i}} < 0 \end{cases} \tag{5}$$

$$S(s) = \frac{n(n-1)(2n+5)}{18}.\tag{6}$$

In Equations (3)–(6), *NDVIi* and *NDVIj* represent the NDVI for years *i* and *j*, respectively, *n* represents the length of the time sequence, and *Sgn* is a symbolic function. The value range of statistic *<sup>Z</sup>* is [−∞, ∞]. At a given significance level α, values of |*Z*|>u1−*α*/2 indicate significant change in the studied sequence at the *α* level. u is a statistic. Here, *Z* = ±2.58 is the significance level of *α* = 0.01, *Z* = ±0.96 is the significance level of *α* = 0.05, and *Z* = ±1.65 is the significance level of *α* = 0.10. In this study, a significance level of *α* = 0.05 was taken to assess the significance of the transformation trend of the NDVI time sequence.

### 2.3.4. Correlation Analysis

The correlation coefficient between the NDVIs and the climate factors based on a pixel was calculated, as follows:

$$R\_{xy} = \frac{\sum\_{i=1}^{n} (\mathbf{x}\_i - \overline{\mathbf{x}})(y\_i - \overline{y})}{\sqrt{\sum\_{i=1}^{n} (\mathbf{x}\_i - \overline{\mathbf{x}})^2} \cdot \sqrt{\sum\_{i=1}^{n} (y\_i - y)^2}},\tag{7}$$

where variable *i* represents the annual serial number, *n* is assigned the value of 21, *xi* represents NDVI data for the *i*-th year, *yi* is the climate factor for the *i*-th year, *x* and *y* are the mean values of variables *x* and *y*, respectively, and *Rxy* is the correlation coefficient between the NDVI and the climate factors. The Mann–Kendall test was used to assess the significance of *Rxy* and to determine its positive and negative correlation thresholds. The analysis was based on the correlation between the annual NDVI, precipitation, and temperature of each pixel.

### **3. Results and Analysis**

### *3.1. Land Use and Land Cover Change Analysis*

As shown in Figure 3, during the study period, snow cover and glaciers were primarily concentrated in high-elevation regions (>3500 m) in the Qilian Mountains. Water bodies included three inland rivers and a number of reservoirs. Forest mainly occupied areas at elevations of 2300–3300 m in the Qilian Mountains. Farmland was mainly distributed in the surroundings of the middle and lower reaches of the rivers. Grassland was divided into natural grassland and artificially planted grassland; the former was located primarily in upstream regions, while the latter was located in middle and lower parts of the river basins. Bare land or desert land covered the largest area in the Hexi Corridor. From 2000 to 2020, the area of land use change was 6.87% of total area.

**Figure 3.** Land use in the Hexi Corridor: (**a**) 2000, (**b**) 2020 and (**c**) changed area.

Table 1 presents a transformation matrix of land use change in the Hexi Corridor from 2000 to 2020, showing the percentage of the area occupied by each land use type. In 2000, the areas of bare land or desert, grassland, farmland, and forest accounted for 68.08%, 21.60%, 5.75%, and 3.05% of the total area, respectively, whereas the area of water bodies, residential land, and snow or glaciers accounted for 0.59%, 0.47%, and 0.46% of the total area, respectively (Table 1). In comparison with 2000, the areas of bare land or desert, grassland, forest, and snow or glaciers in 2020 had decreased to 67.28%, 21.50%, 3.04%, and 0.45% of the total area, respectively; whereas the areas of farmland, residential land, and water bodies had increased to 6.40%, 0.68%, and 0.64% of the total area, respectively. The increases in the area of farmland and residential land represent changes of 0.65% and 0.21%, respectively, whereas the decrease in the area of bare land or desert represents a change of 0.8%. Over the past 21 years, 1.73% of the total area of the Hexi Corridor has been converted from bare land or desert to grassland, while 1.69% of the total area has been converted from grassland to bare land or desert, 0.64% of the total area has been converted from bare land or desert to farmland, and 0.17% of the total area has been converted from farmland to bare land or desert. The interconversion between bare land or desert and forest was equivalent.


**Table 1.** Transfer matrix of land use types in the oases of the Hexi Corridor from 2000 to 2020.

As shown in Table 2, the total NP increased by 15.86, but AREA decreased by 13.69%. Taken together, considering the increase of LSI and SHDI, and the decrease of AI and CON-TAG, it can be concluded that the connectivity of different landscape patches decreased, indicating an obvious landscape fragmentation.

**Table 2.** Landscape pattern index of the entire Hexi Corridor.


*3.2. Characteristics of Vegetation Patterns and Analysis of Vegetation Dynamics in the Hexi Corridor*

Given that the NDVI can reflect the overall condition of regional vegetation, the spatiotemporal distributions of annual maximum NDVIs during 2000–2020 were analyzed (Figure 4). The results ranged from 0 to 0.859, and the average value over the 21-year study period was 0.168. Overall, the annual maximum NDVIs increased with time, and the minimum and maximum values appeared in 2001 and 2019, respectively. In 2010, the annual NDVI increased to the mean state during 2000–2020, indicating improved conditions for vegetation growth in the Hexi Corridor. However, clear regional differences were evident. The highest NDVI values were distributed in the southeast, and the lowest values were found in the northwest. The highest annual average and maximum NDVI values were found in the Qilian Mountains in the southeast of the Hexi Corridor, while the lowest annual average and maximum NDVI values were found in the desert of Jiuquan, which might reflect climatic influences, especially that of precipitation.

**Figure 4.** Spatiotemporal distribution of growing season NDVI during 2000–2020.

The land use of regions with higher annual average and maximum NDVI values mainly comprised forest, grassland, and farmland, whereas that of the areas with lower NDVI values was mainly desert vegetation. Over the 21-year study period, the average NDVI of forest, grassland, farmland, bare land, and desert was 0.618, 0.486, 0.515, and 0.112, respectively, indicating that the environment for vegetation growth was better in forest, farmland, and grassland areas than in desert areas.

### *3.3. Quantitative Analysis of Driving Factors of Interannual Variation of NDVI in the Hexi Corridor*

### 3.3.1. Spatial Patterns of Climatic Effects on Vegetation Changes

The meteorological data were gridded using multivariate regression analysis and Kriging interpolation [24], then used to calculate the correlation with NDVI.

As shown in Figure 5a, except for bare land or desert areas, there was positive correlation between the NDVI and precipitation. Specifically, 41.36% of oasis vegetation areas were correlated positively (*p* < 0.05), of which 13.02% was significant at the 0.01 level. The areas of positive correlation included grassland (79.21%), farmland (12.82%), and forest (4.06%), distributed mainly in the south of Wuwei and Jiuquan counties, to the west of Jinchang, and in the Qilian Mountains. A potential reason for this is that the increasing growth of farmland and grassland vegetation in this area, together with the amount of precipitation, had a substantial impact on the regional NDVI. Only 0.03% of areas had a negative correlation.

**Figure 5.** Correlation analysis of NDVI and (**a**) precipitation and (**b**) temperature changes in the Hexi Corridor during 2000–2020.

In terms of the temperature effect (Figure 5b), 5.38% of oasis vegetation areas were correlated negatively (*p* < 0.05), of which 1.65% was significant at the 0.01 level. The areas of negative correlation were bare land and desert (53.23%), grassland (29.2%), and farmland (16.76%); mainly distributed in the northwest of Minqin, the southeast of Zhangye, and the southwest of Jiuquan. This demonstrated that the vegetation change over the entire region was affected more by precipitation than by temperature.

### 3.3.2. Effects of Topographic Factors on Spatial Patterns of Vegetation Changes

The changes of vegetation at the different elevations during 2000–2020 were discussed based on DEM data. Overall, the NDVI presented an upward trend, but with obvious spatial heterogeneity (Figure 6). For 5.66% of the total area, the NDVI increased significantly (*p* < 0.05), of which 1.85% was significant at the 0.01 level. For 1.1% of the total area, the NDVI decreased significantly (*p* < 0.05), of which 0.21% was significant at the 0.01 level. Generally, the NDVI increased at the edges of the Hexi Corridor oases, while it decreased in southern parts of the city of Dunhuang, central parts of the city of Zhangye, and central parts of Wuwei and Minqin counties. As for the analysis of land use type, 55.65% and 33.79% of the increasing NDVI areas were farmland and planted grassland, while 73.05%

and 13.84% of the decreasing NDVI areas were farmland and building land. The changes to NDVI were in flat areas with >3000 m and <2000 m, but fluctuated greatly in the anthropic area from 2000 to 3000 m. Overall, this was mainly influenced by human activities.

**Figure 6.** Trend of the NDVI in the Hexi Corridor during 2000–2020.
