## 3.1 Vegetation NDVI classification

According to the equidistant method, the annual NDVI data were divided into five vegetation coverage grades: low vegetation cover (L-NDVI) (NDVI ≤ 0.2), medium-low vegetation cover (ML-NDVI) (0.2 < NDVI ≤ 0.4), medium vegetation cover (M-NDVI) (0.4 < NDVI < 0.6), medium-high vegetation cover (MH-NDVI) (0.6 < NDVI ≤ 0.8) and high vegetation cover (H-NDVI) (NDVI > 0.8) according to the value of NDVI, so as to better analyze the dynamic changes of NDVI (Shi et al. 2023).

## 3.2 Coefficient of variation

The coefficient of variation (CV) is often used to indicate the degree of fluctuation of time series data, which can indicate the vulnerability of regional ecosystem. The formula as follows:

$$\text{C}\text{V}=\frac{1}{\stackrel{-}{x}}\sqrt{\frac{1}{n-1}\sum _{i=1}^{n}{({x}_{i}-\stackrel{-}{x})}^{2}}$$

1

Where, CV is the coefficient of variation, *i* is time series, *x**i* is the pixel value of the *i*th year, and \(\stackrel{-}{x}\) is the mean value of the pixel values of all years in the study period.

Taking the interannual NDVI as the input value, the larger CV values indicates that the vegetation coverage changes drastically and the ecological environment is fragile. The smaller CV values indicates that the fluctuation of vegetation coverage is small and the anti-interference ability of ecological environment is strong (Han et al.2022). The coefficient of variation is classified, as shown in Table 1.

Table 1

Classification of CV | Value of CV |

High fluctuation change | CV > 0.3 |

Relatively high fluctuation change | 0.2<CV ≤ 0.3 |

Moderate fluctuation change | 0.15<CV ≤ 0.2 |

Relatively low fluctuation change | 0.10<CV ≤ 0.15 |

Low fluctuation change | CV ≤ 0.1 |

## 3.3 Theil-Sen Median and Mann-Kendall test

Theil-Sen Median (Sen) and Mann-Kendall test (MK) are widely used in meteorological and hydrological research (Dagnachew et al. 2020; Li et al. 2021). The Sen analysis can effectively avoid outliers or measurement errors, and is not interfered by outliers. At the same time, it is more suitable than the linear regression method to study the trend of vegetation change. The calculation formula is as following:

$$\rho =median\left\{\frac{{x}_{j}-{x}_{i}}{j-i}\right\},1<i<j<n$$

2

Where, *ρ* is the slope of time series data, when *ρ* > 0, NDVI showed an increasing trend, otherwise, it shows a downward trend. Mann-Kendall test was used to determine the significance of NDVI trend.

Mann-Kendall test is a non-parametric test method. It does not require the sample to follow a specific distribution. It is less disturbed by outliers and is more suitable for sequential variables (Fan et al. 2020). The Mann-Kendall test was used primarily to test the significance through the Z value. In this study, the significance of NDVI change trend at 0.05 confidence level was selected.

When |Z| > Z1−α/2, NDVI shows a significant change trend. Z1−α/2 is the corresponding value of the standard normal distribution function at the α level. The trend of NDVI is categorized into six levels based on *ρ* and |Z|, namely very significantly elevated(*ρ* ≥ 0, |Z| > 2.58, P < 0.05), Very significant reduction(*ρ* < 0, |Z| > 2.58, P < 0.05), Significantly higher(*ρ* ≥ 0, |Z| > 1.96, P < 0.05), Significantly lower(*ρ* < 0, |Z| > 1.96, P < 0.05), No significant elevation(*ρ* ≥ 0, |Z| > 1.96, P > 0.05), No significant reduction(*ρ* < 0, |Z| ≤ 1.96, P > 0.05).

## 3.4 Hurst Index

The Hurst index is able to quantitatively characterize the long-term dependence of time-series information and forecast its future evolution (Zhong et al. 2022). The principle is as follows:

$${\stackrel{-}{NDVI}}_{\left(\tau \right)}=\frac{1}{\tau }{\sum }_{t=1}^{\tau }{NDVI}_{\left(t\right)}，{\tau }=\text{1,2},\dots ,\text{n}$$

3

\({X}_{\left(t,\tau \right)}={\sum }_{t=1}^{t}\left({NDVI}_{\left(t\right)}-{\stackrel{-}{NDVI}}_{\left(\tau \right)}\right)\) , 1\(\le t\le \tau\) (4)

$${R}_{\left(\tau \right)}=\begin{array}{c}max\\ 1\le t\le \tau \end{array}{X}_{\left(t,\tau \right)}-\begin{array}{c}min\\ 1\le t\le \tau \end{array}{X}_{\left(t,\tau \right)}，{\tau }=\text{1,2},\dots ,\text{n}$$

5

$${S}_{\left(\tau \right)}={\left[\frac{1}{\tau }{\sum }_{t=1}^{\tau }{\left({NDVI}_{\left(t\right)}-{\stackrel{-}{NDVI}}_{\left(t\right)}\right)}^{2}\right]}^{\frac{1}{2}},{\tau }=\text{1,2},\dots ,\text{n}$$

6

$$\frac{{R}_{\left(\tau \right)}}{{S}_{\left(\tau \right)}}={\left(c\tau \right)}^{H}$$

7

Where, (3) to (7) are the series of NDVI mean, cumulative deviation, extreme deviation and standard deviation, respectively. H is the Hurst index, which can be found by linear regression analysis.

When 0 < H < 0.35, it shows strong inverse continuity, when 0.35 < H < 0.5, it shows weak inverse continuity, when 0.5 < H < 0.65, it shows weak positive continuity, when 0.65 < H < 1, it shows strong positive continuity.

## 3.5 GeoDetector

GeoDetector is a spatial statistical method for detecting the spatial heterogeneity and drivers of geographic objects, and the model consists of four modules, namely Factor_detector, Risk_detector, Interaction_detector, and Ecological_detector. The Factor_detector quantitatively describes the relative importance of the influencing factors, and it measures the ability of the independent variable to explain the dependent variable by constructing a *q* statistic. Combining the geographical characteristics of the study area and considering the natural as well as anthropogenic factors the selected influencing factors are given in Table 2 (Zhu et al. 2020).

Table 2

Indicators of detection factors

Impact Factor | Indicators | Unit Symbol |

X1 | Temperature | ℃ |

X2 | Precipitation | mm |

X3 | Annual sunshine hours | h |

X4 | Relative humidity | %RH |

X5 | Evaporation | mm |

X6 | DEM | m |

X7 | Slope | ° |

X8 | Slope direction | - |

X9 | Landform type | - |

X10 | Population density | person/km2 |

X11 | Ground temperature | ℃ |

X12 | Wind speed | m/s |

X13 | Land use type | - |

GeoDetector include four sub-detectors:

(1) Factor_detector

It is utilized to detect spatially stratified heterogeneity in the dependent variable and also to reflect the extent to which the respective variable explains the spatially stratified heterogeneity in the dependent variable. The larger the value of *q*, The formula is as follows:

$$q=1-\frac{\sum _{h=1}^{L}{N}_{h}{\sigma }_{h}^{2}}{N{\sigma }^{2}}$$

8

where, *q* takes the value of 0 ~ 1; *N**h* is the number of sample units in layer *h*, *N* is the number of sample units in the whole district, *σ*h2 is the variance of layer *h* with respect to the value of the dependent variable Y, *σ*2 is the variance of layer *h* with respect to the value of the dependent variable Y. Combining the existing studies, the explanatory power was categorized into four categories, namely strong influence (0.5 ≤ *q* ≤ 1), moderate influence (0.3 ≤ *q* < 0.5), weak influence(0.1 ≤ *q* < 0.3), and none (0 ≤ *q* < 0.1) (Chen et al. 2022).

(2) Interaction_detector

Interaction_detector is used to determine the interaction of different risk factors on dependent variables. The type of interaction between the two factors is given in Table 3.

Table 3

interaction | interactive relationship |

*q*(X1∩X2) < min[*q*(X1), *q*(X2)] | Nonlinear weakening |

min[*q*(X1), *q*(X2)] < *q*(X1∩X2) <max[*q*(X1),*q*(X2)] | Single-factor nonlinear attenuation |

*q*(X1∩X2) > max[*q*(X1), *q*(X2)] | Two-factor enhancement |

*q*(X1∩X2) = *q*(X1) + *q*(X2) | Dependency |

*q*(X1∩X2) > *q*(X1) + *q*(X2) | Nonlinear enhancement |

(3) Risk_detector.

To determine whether there is a significant difference in the mean value of the attributes between the two sub-regions, the *t* statistic is used for measurement:

$${t}_{{y}_{h=1}-{y}_{h=2}}=\frac{{\stackrel{-}{Y}}_{h=1}-{\stackrel{-}{Y}}_{h=2}}{{\left[\frac{var\left({\stackrel{-}{Y}}_{h=1}\right)}{{n}_{h=1}}+\frac{var\left({\stackrel{-}{Y}}_{h=2}\right)}{{n}_{h=2}}\right]}^{1/2}}$$

9

where, \({\stackrel{-}{Y}}_{h}\) is the mean N DVI attribute value of the region; *n**h* is the number of samples in the subregion.

(4) Ecological_detector

The *F* statistic is used to refer to the variability of the effect of the two factors in the independent variable on the spatial distribution of the attributes of the dependent variable.

$$F=\frac{{N}_{{X}_{1}}\left({N}_{{X}_{2}}-1\right){A}_{{X}_{1}}}{{N}_{{X}_{2}}\left({N}_{1}-1\right){A}_{{X}_{2}}}$$

10

$${A}_{{X}_{1}}=\sum _{h=1}^{{L}_{1}}{N}_{h{\sigma }_{h}^{2}}$$

11

$${A}_{{X}_{2}}=\sum _{h=1}^{{L}_{2}}{N}_{h}{\sigma }_{h}^{2}$$

12

Where, \({N}_{{X}_{1}}\)and \({N}_{{X}_{2}}\) are the sample sizes of the dependent variables *X**1* and *X**2*, respectively, \({A}_{{X}_{1}}\) and \({A}_{{X}_{2}}\) are the sums of the within-stratum variances of the strata formed by *X**1* and *X**2*, respectively, *L**1*, *L**2* are the number of strata of the independent variables, respectively.