Topographic heterogeneity regulates nonlinear vegetation responses to climate and human activities in China’s ecological transition zone

1 Introduction

Vegetation is a fundamental component of terrestrial ecosystems. It regulates the carbon and water cycles, facilitating energy exchange between the Earth’s surface and the atmosphere. This process connects the atmosphere, hydrosphere, and lithosphere, underscoring vegetation’s crucial role in sustaining the functions of terrestrial ecosystems (Piao et al., 2020; Han et al., 2024; Chen et al., 2025; Shang et al., 2026). Furthermore, vegetation responds rapidly to environmental changes, making shifts in vegetation reliable indicators of ecosystem dynamics in the context of changing climatic conditions (Kong et al., 2017; Zhang X. Y. et al., 2024; Sun et al., 2026). Additionally, alterations in vegetation have direct implications for the global carbon cycle, biodiversity, and ecosystem services (Hautier et al., 2015; Spohn et al., 2023).

Satellite remote sensing has significantly improved the monitoring of vegetation dynamics on both global and regional scales by providing continuous observations over extensive spatial and temporal ranges. Among the various indices available, the normalized difference vegetation index (NDVI) is one of the most widely utilized indicators for tracking vegetation activity. Its effectiveness in characterizing vegetation cover and growth conditions contributes to its popularity (Huang et al., 2021; Yang et al., 2024; Zhao and Qu, 2024). Consequently, long time-series NDVI products derived from sensors such as Landsat, AVHRR, MODIS, and SPOT-VGT serve as a crucial data foundation for studying vegetation responses to environmental changes and human disturbances (Albarakat and Lakshmi, 2019). In the context of global warming and increasing human activities, research focused on vegetation changes and their driving factors at both global and regional scales, based on long-term NDVI datasets, has emerged as a key area of interest. This research holds significant scientific and practical importance for understanding ecosystem stability, conserving biodiversity, assessing carbon sinks, and developing adaptation strategies at both global and regional levels (Dutta and Markonis, 2024).

Numerous studies have demonstrated that global vegetation changes exhibit significant spatial variability over the past few decades (Zhu et al., 2016; Zeng et al., 2018; Piao et al., 2020; Cai et al., 2025). For example, vegetation has markedly improved in the mid- to high latitudes of the Northern Hemisphere. This improvement is largely attributed to a lengthening growing season resulting from warming, which has fostered increased plant growth. China and India have been pivotal in this greening trend, primarily due to extensive ecological restoration initiatives and effective land management practices (Chen et al., 2019).

Vegetation changes in China exhibit significant spatial variability, primarily due to diverse climatic conditions, variations in water and heat availability, and the presence of different ecosystem types across regions (Yao et al., 2019; Zhang and Ye, 2020). Climatic factors, including air temperature, precipitation, and vapor pressure deficit (VPD), are crucial in regulating vegetation productivity. These factors influence photosynthesis, phenology, and water availability (Yuan et al., 2019; Wu et al., 2021; Zhang M. et al., 2024). Additionally, human activities—such as land-use changes, urban expansion, and ecological restoration—have emerged as significant drivers of vegetation dynamics (Li et al., 2020; Tian et al., 2021).

The Qinling–Daba Mountains represent one of central China’s most significant ecological transition zones, acting as a natural boundary between northern and southern China in terms of climate, ecology, and biogeography. This region is influenced by both the East Asian monsoon and the climate system of the Qinghai-Tibet Plateau. It serves as a crucial ecological security barrier that connects the Yangtze and Yellow River basins, while also providing essential ecosystem services, including biodiversity preservation and water conservation (Zhang Y. et al., 2022). In mountainous areas, complex topography significantly affects the nonlinear responses of vegetation by altering mountain climate, hydrological processes, and the intensity of human activities. Consequently, the effects of climatic factors may be obscured by topography-driven resource distribution or the extent of human disturbance (Sanz et al., 2021). In ecologically transitional zones marked by intricate topography, the relative contributions of various drivers to vegetation change remain a topic of debate (Wang et al., 2023). Amid global change, climatic conditions in the Qinling–Daba Mountains ecological transition zone have undergone significant shifts, exhibiting north–south differentiation influenced by topography. Simultaneously, human activities have evolved within the constraints imposed by the terrain (Zhang X. et al., 2022). In this context, changes in vegetation within the region are influenced by topography and are particularly sensitive to climate change and human disturbances. However, previous studies have largely neglected the moderating effects of topographic heterogeneity—such as altitude, slope, and aspect—on vegetation responses (Xu S. et al., 2020; Tian and Tian, 2026).

This study investigates the Qinling–Daba Mountains, a key ecological transition zone, by integrating long-term (2000–2022) 30-meter Landsat-derived NDVI data with various datasets on climate, topography, and human activities. The aim is to explore vegetation dynamics and their driving factors. Specifically, we utilize Theil–Sen estimates, Mann–Kendall tests, and the Hurst index, in conjunction with a Bayesian-optimized explainable machine learning approach (XGBoost–SHAP), to: (1) characterize the spatial and temporal patterns of vegetation changes and their persistence; (2) assess the influence of topographic heterogeneity on vegetation responses to climatic and anthropogenic drivers; and (3) identify nonlinear threshold responses and interaction mechanisms among multiple environmental factors. This study goes beyond merely ranking the importance of environmental variables; it seeks to elucidate how and under what environmental conditions vegetation responds to climate and human-induced changes. In doing so, it provides a transferable analytical framework for understanding vegetation–environment interactions and supports ecosystem restoration efforts in ecologically sensitive mountainous regions.

2 Overview of study area

The Qinling–Daba Mountain Region (30°30′–34°37′N, 103°44′–113°13′E) is situated at the intersection of Shaanxi, Gansu, Sichuan, Chongqing, Hubei, and Henan in central China. This region encompasses the Qinling Mountains, the Han River Valley, and the Daba Mountains, serving as a significant ecological and climatic transition zone between northern and southern China. Spanning approximately 304,986 km², it is characterized by complex mountainous terrain and notable elevation gradients (Figure 1). The region displays pronounced latitudinal and altitudinal climatic gradients. The northern areas experience a warm-temperate semi-humid to semi-arid monsoon climate, while the southern areas are characterized by a subtropical monsoon climate. Mean annual temperatures range from 12 to 16 °C, and annual precipitation varies from 709 to 1,500 mm (Lian et al., 2024). These climatic and topographic gradients contribute to significant vegetation heterogeneity, which includes temperate deciduous broad-leaved forests in mountainous regions and mixed evergreen-deciduous broad-leaved forests in river valleys. This combination of pronounced climatic, topographic, and vegetation diversity makes the Qinling–Daba Mountains an ideal natural laboratory for examining how topography influences vegetation responses to climatic and anthropogenic factors (Zhang et al., 2021; Yao, 2024; Zou et al., 2026).

3 Data and methods

3.1 Data and its pre-processing

3.1.1 NDVI dataset

To characterize long-term vegetation dynamics, we utilized annual maximum NDVI at a spatial resolution of 30 m for the period from 2000 to 2022. The dataset was generated by the Remote Sensing Team for Land Use and Global Change at the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences (Yang et al., 2019). Imagery from Landsat 5, 7, and 8 was processed using the Google Earth Engine platform, where observations contaminated by clouds and shadows were removed prior to NDVI reconstruction. Subsequently, temporal interpolation and Savitzky–Golay filtering were applied to smooth the time series and produce continuous annual maximum NDVI products for China. For the subsequent driver analysis, NDVI was aggregated from the original 30 m resolution to a common 1-km analysis grid using mean aggregation. This step ensured spatial consistency with the climate and socioeconomic predictors.

3.1.2 Climate and related data

Climate-related variables included air temperature (TMP), precipitation (PRE), and potential evapotranspiration (PET). All climatic datasets were clipped to the study area and converted into annual time-series data. Annual mean or cumulative values were then calculated to represent regional climatic background conditions. These climatic variables primarily served to evaluate the impact of climate change on vegetation dynamics and were utilized as explanatory factors in driver analyses. To further assess vegetation responses to water stress, vapor pressure deficit (VPD) and soil moisture (SM) were incorporated into the analysis. Following Liu Y. Y. et al. (2024), VPD was derived from the ERA5-Land reanalysis dataset using the 2 m air temperature and 2 m dew point temperature variables. VPD characterized atmospheric drought conditions, while soil moisture indicated surface water availability. These indicators were employed to examine the effects of water stress on vegetation variability. Ultimately, all the aforementioned factors (TMP, PRE, PET, VPD, and SM) were employed as driving variables in an XGBoost-SHAP analysis. Prior to this analysis, all raster layers were resampled to a uniform spatial resolution of 1 km using bilinear interpolation and clipped to the study area boundary.

3.1.3 Other datasets

Topographic information was derived from the Shuttle Radar Topography Mission Digital Elevation Model (SRTMGL1), released by the U.S. Geological Survey, with a spatial resolution of 30 m. Utilizing the DEM data, we extracted terrain-related factors, such as slope and aspect, to assess the impact of topographic conditions on vegetation distribution and change patterns. For the driver analysis, continuous topographic variables were aggregated to a uniform 1-km analysis grid. Elevation and slope were calculated using area-weighted mean values derived from fine-resolution pixels within each 1-km grid cell. To accurately represent aspect, we employed circular vector-weighted averaging, which accommodates its circular angular nature (0–360°) and mitigates biases associated with simple arithmetic averaging. Human activity indicators mainly included population density and land use type. Population density data were obtained from the LandScan dataset developed by the Oak Ridge National Laboratory of the U.S. Department of Energy, with a spatial resolution of 1 km. Before conducting the driver analysis, these values were transformed using a logarithmic function with a constant offset [lg(Pop+1) + 1]. This transformation aimed to reduce data skewness and ensure numerical stability when these values were incorporated as predictors in the XGBoost model. Land use data, taken from the 30-m resolution dataset of Yang and Huang (2026), were aggregated to the 1-km grid by majority rule and incorporated as a driver in the XGBoost-SHAP model. These variables were incorporated to explore the effects of anthropogenic activities and their interactions with environmental factors on vegetation dynamics. Additionally, the study employed a 1:1,000,000 vegetation type dataset for China to investigate vegetation change characteristics across different ecosystem types. Detailed original data sources are shown in Table 1.

Table 1. Data sources and detailed descriptions.

3.2 Methods

The framework of this study is shown in Figure 2. We focused on the Qinling–Daba Mountains, a representative climate transition zone and ecological barrier in China, as our study area. We integrated multi-source datasets to examine vegetation dynamics and their driving mechanisms from 2000 to 2022. The analytical framework comprised three main stages: (1) assessing long-term vegetation dynamics and persistence; (2) correcting for terrain-related spatial biases; and (3) quantifying the effects and nonlinear responses of climatic, topographic, and anthropogenic factors using Bayesian-optimized XGBoost-SHAP.

3.2.1 Methods for Analyzing spatiotemporal trends in vegetation

3.2.1.1 Theil–Sen median trend analysis

The Theil–Sen approach (Sen slope estimate) estimates the median slope among all pairwise combinations within a time series and is less sensitive to abnormal values or measurement uncertainties, making it suitable for long-term vegetation monitoring (Zhang S. et al., 2022). The corresponding calculation is presented in Equation 1.

S NDVI = Median ( NDVI j − NDVI i j − i ) , ( 2000 ≤ i < j ≤ 2022 ) (1)

In the equation, Median() denotes the midpoint; NDVI j and NDVI i denote NDVI data for different years. When S NDVI>0, it indicates an upward trend in NDVI; conversely, it indicates a downward trend in NDVI.

3.2.1.2 Mann–Kendall test

To evaluate the statistical significance of vegetation trends, the nonparametric Mann–Kendall test was further applied. Since this method is independent of data distribution assumptions and exhibits strong resistance to outliers, it has been adopted in ecological and environmental time-series analyses (Al Mashoudi et al., 2025).

Setup { NDVI i } , i = 2000 , 2001 , … , 2022

Define Z as the statistic, whose formula is given in Equation 2.

Z = { S − 1 s ( S ) , S > 0 0 , S = 0 S + 1 s ( S ) , S < 0 (2)

In the equation, S = ∑ j = 1 n − 1 ∑ i = j + 1 n sgn ( NDVI j − NDVI i ),s ( S ) = n ( n − 1 ) ( 2 n + 5 ) 18, sgn is a symbolic function; its formula is given in Equation (3).

sgn ( NDVI j − NDVI i ) = { 1 , NDVI j − NDVI i > 0 0 , NDVI j − NDVI i = 0 − 1 , NDVI j − NDVI i < 0 (3)

In the equation, NDVI j and NDVI i denote the NDVI values of pixel i in years i and j, respectively, and n denotes the length of the time series.

The range of the statistic Z is (− ∞ , + ∞). At a given significance level, when ∣ Z ∣ > u 1 − α / 2, it indicates that there is a significant change in the study series at α level. Generally, α is taken as 0.05; in this paper, we assess the significance of the trend in the NDVI time series at the α = 0.05 confidence level.

3.2.1.3 Hurst exponent analysis

The Hurst index analysis is an effective method for quantitatively describing long-range dependence in time series and has been widely applied in the field of vegetation change (Sanz et al., 2021; Wang et al., 2023; Zhang et al., 2025). Its basic principle is as follows: for a time series { NDVI ( t ) } . t = 1 , 2 , … , n, a mean sequence is defined, the calculation of which is shown in Equation 4.

NDVI ¯ ( T ) = 1 T ∑ t = 1 T NDVI ( T ) (4)

The cumulative deviation (X ( t , T )), range (R ( T )), and standard deviation (S ( T )) of this time-series data can be expressed using Equations 5–7, respectively.

X ( t , T ) = ∑ t = 1 T ( NDVI ( t ) − NDVI ¯ ( T ) ) (5)

R ( T ) = max 1 ≤ t ≤ T X ( t , T ) − min 1 ≤ t ≤ T X ( t , T ) (6)

S ( T ) = [ 1 T ∑ t = 1 T ( NDVI ( t ) − NDVI ( T ) ) 2 ] 1 2 (7)

In particular, 1 ≤ t ≤ T, T = 1 , 2 , … , n.

For the ratio R ( T ) S ( T ) ≜ R S, if the following relationship holds: R S ∝ T H, this indicates that the analyzed time series data exhibits the Hurst phenomenon, and H is referred to as the Hurst index. The value of H can be obtained by fitting Equation 8 using the least squares method.

log ( R S ) n = a + H × log ( n ) (8)

The Hurst exponent was used to evaluate whether the NDVI time series exhibits persistence in its temporal evolution. Specifically, H > 0.5 indicates persistent behavior, suggesting that future changes are likely to continue in the same direction as historical trends. A larger H value corresponds to stronger persistence. When H = 0.5, the series behaves as a random process without long-term dependence. In contrast, H < 0.5 indicates anti-persistent behavior, implying that future variations tend to reverse previous trends, with smaller H values representing stronger anti-persistence. This study predicts the persistence of future vegetation changes in the Qinling–Daba Mountain region by combining Hurst index analysis with Theil–Sen median trend analysis.

3.2.2 Terrain area difference correction

Topographic heterogeneity may introduce bias when comparing vegetation change characteristics across different terrain conditions. To reduce the influence of unequal area distributions among topographic classes, a terrain area correction coefficient (K) was introduced in this study. Based on the proportional changes in vegetation improvement, degradation, and stable areas under different elevation, slope, and aspect conditions, the coefficient was used to characterize the relationship between topography and vegetation dynamics (Li et al., 2023). The calculation of K is given in Equation 9.

K = S ie S e × ( S i S ) − 1 (9)

Here, S ie represents the area of a specific vegetation change trend type i (vegetation increase, vegetation decrease, or vegetation stability) under specific topographic conditions e (a given range of elevation, slope, or aspect); S e represents the total area of topographic condition e; S ie S e represents the proportion of area occupied by vegetation change trend type i under topographic condition e; S i represents the total area of vegetation change trend type i; and S represents the total area of the study region.

When K > 1, it indicates that a specific vegetation change trend type is dominant under e topographic conditions; K = 1 indicates that a specific vegetation change trend type is evenly distributed under e topographic conditions; K < 1 indicates that a specific vegetation change trend type is not dominant under e topographic conditions.

3.2.3 Bayesian optimization-based XGBoost-SHAP

XGBoost is a machine learning algorithm, which employs the concept of gradient boosting decision trees to improve predictive performance by iteratively training decision tree models. Compared to traditional machine learning algorithms, this algorithm enhances the model’s ability to identify optimal solutions (Chen and Guestrin, 2016). The algorithm can be expressed by Equation 10.

y i ̂ = ∑ k = 1 k f k ( x i ) , f k ∈ F ( i = 1 , 2 , … , n ) (10)

Here, F = { f ( x ) = ωq ( x ) } ( q : Rm → { 1 , 2 , … , T } , ω ∈ RT ) denotes the decision trees in the ensemble model of classification and regression trees; q maps samples to the leaf nodes of the tree structure; T is the number of leaf nodes; and ω is the weight of the leaf nodes.

The objective function in the XGBoost algorithm consists of two parts: an error function and a complexity function, as shown in Equation 11.

Obj = ∑ i = 1 n ( y i − y i ̂ ) 2 + ( γT + 1 2 γ ∑ j = 1 T ω j 2 ) (11)

Here, 𝛾𝑇 represents the L1 regularization term, and 1 2 γ ∑ j = 1 T ω j 2 represents the L2 regularization term. Building on this, the predicted values from the tth-round model are incorporated, and a new tth tree is introduced. The objective function is approximated using a second-order Taylor expansion, ultimately yielding the optimized objective function.

Bayesian optimization was adopted to determine the optimal hyperparameter combinations for the machine learning models. Compared with traditional parameter search approaches, Bayesian optimization can iteratively update the search strategy according to previous evaluation results, thereby improving optimization efficiency (Shahriari et al., 2015; Alibrahim and Ludwig, 2021). During the optimization process, the algorithm continuously explored the hyperparameter space and identified the parameter combination that achieved the best model performance (Sun et al., 2021).

SHAP is a method of interpreting the predictions of machine learning models that is based on the Shapley values in game theory (Lundberg and Lee, 2017). It is often used to quantify the contribution of each feature to the model’s predictions and to distinguish between positive and negative effects. It can overcome the inherent “black-box” nature of the XGBoost model, thereby improving its interpretability and transparency.

In this study, the NDVI was selected as the dependent variable. The independent variables included TMP, PRE, PET, SM, VPD, elevation (DEM), slope, aspect, population density (Pop), and land use type (LULC). To analyze these variables, three independent sampling datasets were created using the Subset Feature tool in ArcGIS. Initially, 2,000 points were randomly selected from the entire study area to represent overall vegetation conditions. Subsequently, 1,000 points were chosen from regions exhibiting significant increases in NDVI, while another 1,000 points were selected from areas demonstrating significant decreases, as identified by the Sen–MK analysis. The Multivalue to Point tool was then employed to construct a time-series dataset. The XGBoost regression model, combined with SHAP interpretation, was employed to quantify the effects of multiple factors on NDVI variation. To prevent spatial leakage, a 10-km spatial-block grouping strategy was utilized for dataset partitioning. Specifically, blocks, rather than individual sampling points, were randomly divided into training (80%) and test (20%) subsets. All temporal records associated with the same Block_ID and POINT_ID were retained entirely within either the training or test set, ensuring spatial independence between the two datasets. Hyper-parameter tuning was conducted using Bayesian optimization within the training set, accompanied by spatially-constrained Group-KFold cross-validation based on Block_ID to address spatial autocorrelation. Following the identification of optimal hyper-parameters, the final XGBoost model was refitted using the complete training dataset. Model performance was assessed independently on the held-out test dataset through the coefficient of determination (R²), root-mean-square error (RMSE), and mean absolute error (MAE). This modeling workflow was executed separately for the overall samples, significantly increasing NDVI samples, and significantly decreasing NDVI samples. Additionally, SHAP analysis was applied to interpret variable importance and the marginal response relationships between predictors and NDVI. SHAP importance rankings, honeycomb plots, and factor dependency plots were utilized for interpretable analyses of the drivers of vegetation change.

4 Results

4.1 Spatiotemporal characteristics of NDVI in the Qinling–Daba Mountains during 2000–2022

4.1.1 Temporal trends in vegetation

Annual NDVI in the Qinling–Daba Mountains exhibits a significant upward trend over the past 23 years (Figure 3), with a linear increase rate of approximately 0.0051 yr.⁻¹, indicating continuous vegetation greening. During this period, NDVI fluctuates between 0.72 and 0.76 from 2000 to 2010, followed by a decline to a minimum value of approximately 0.70 in 2012. Thereafter, NDVI increases rapidly and enters a phase of accelerated greening, reaching 0.84 in 2022, the highest value recorded during the study period.

4.1.2 Spatial distribution and trends of vegetation

Vegetation cover in the Qinling–Daba Mountains remains at a relatively high level and exhibits significant spatial heterogeneity (Figure 4). Areas with high vegetation cover account for 77.15% of the total area and are mainly distributed in contiguous patches along the main mountainous regions and southern parts of the Qinling–Daba Mountains. Areas with medium-low and low vegetation cover together account for only 1.11% of the total area and are scattered in fragmented patches across the northwestern and northern regions of the Qinling–Daba Mountains.

From 2000 to 2022, vegetation cover exhibited a clear greening trend. The proportion of high vegetation cover areas increases significantly to 90.13% (Table 2). This expansion is mainly driven by the upward transition of medium-high vegetation cover areas since 2000, with a cumulative transfer area of 48,758.68 km². Meanwhile, the proportions of all other vegetation cover classes decrease to varying degrees, indicating widespread and persistent vegetation improvement across the Qinling–Daba Mountains.

Table 2. Vegetation transfer matrix from 2000 to 2022 (unit: km², the numbers in parentheses represent the percentage of the total area.)

4.1.3 Spatiotemporal characteristics of vegetation change trends

Vegetation trends in the Qinling–Daba Mountains during 2000–2022 were dominated by widespread greening with pronounced spatial heterogeneity. Areas showing significant and extremely significant increases formed extensive contiguous patches across the major mountainous regions (Figure 5a), whereas vegetation degradation occupied less than 8% of the study area and mainly occurred in fragmented patches within low-elevation valleys and urbanized regions. Statistical results further showed that approximately 90.03% of the study area experienced increasing vegetation trends (Figure 5b), including 21.84, 22.35, and 45.84% exhibiting extremely significant, significant, and non-significant increases, respectively. Overall, vegetation dynamics in the Qinling–Daba Mountains were characterized by extensive and persistent greening during the past two decades.

Vegetation dynamics in the Qinling–Daba Mountains were dominated by increasing NDVI trends during 2000–2022. The Theil–Sen slope ranged from −0.0576 to 0.0471, with positive trends occurring across most areas of the region. Hurst analysis further indicated strong persistence in vegetation dynamics (Figures 6a,b), with persistent increase areas accounting for approximately 89.83% of the study area. These areas were mainly distributed across high-altitude mountainous regions, whereas persistent degradation was concentrated in river valleys and urban centers. Anti-persistent areas accounted for less than 0.35% of the region, indicating relatively limited instability in vegetation evolution.

4.2 Spatiotemporal characteristics of main climate variables in the Qinling–Daba Mountains during 2000–2022

4.2.1 Temporal trends in temperature and precipitation

As shown in Figure 7, the climate conditions in the Qinling–Daba Mountains have demonstrated distinct warming and moistening trends over the past 23 years. The annual mean temperature has steadily increased, ranging from 10.4 to 11.4 °C, with an average rate of approximately 0.02 °C yr.⁻¹. Similarly, annual precipitation exhibits a fluctuating upward trend, varying between 650 and 1,050 mm, with an average increase of about 3.2 mm yr.⁻¹. These concurrent increases in temperature and precipitation are likely linked to the sustained greening of vegetation observed in the study area.

4.2.2 Spatiotemporal characteristics of main climate variables trends

The spatial differentiation of temperature and precipitation in the Qinling–Daba Mountains is pronounced. The multi-year average cumulative precipitation reveals a pattern characterized by higher values in the southeast and lower values in the northwest. The southeastern region receives substantial precipitation with a negligible declining trend in interannual variability. In contrast, the northwestern area exhibits a significant increasing trend (Figure 8a). Precipitation levels range from 386.33 to 1,116.46 mm, with a contiguous high-value center in the south that decreases toward the north, resulting in a low-value area in the northern region. This overall pattern illustrates a “high in the south, low in the north” distribution. Topographic factors significantly influence precipitation distribution. The interannual change rate fluctuates between −8.61 and 7.82 mm·yr.⁻¹, underscoring considerable spatial heterogeneity (Figure 8b). In the western region, increasing trends are predominant, with some areas exhibiting extremely significant or significant levels. Conversely, the eastern region shows decreasing trends, although most do not meet significance thresholds. The distribution of areas with significant changes is patchy, revealing a stark spatial contrast between the contiguous zone of extremely significant increases in the west and the decreasing areas in the east. The spatial pattern of multi-year average temperatures displays a gradient of “high in the east, low in the west,” ranging from −9.88 to 18.21 °C, indicating significant spatial differences (Figure 8c). High-temperature areas are primarily concentrated in the southeastern region and the southern river valleys, while low-temperature areas are located in the northwestern high-elevation mountains. The distribution of low temperatures closely corresponds to topographic features, generally decreasing with increasing elevation. The rate of temperature change varies from −0.12 to 0.45 °C·yr.⁻¹, with warming trends prevailing across the entire region; however, the spatial pattern remains notably fragmented (Figure 8d). Significant warming patches are scattered throughout the southeastern and northwestern regions, while isolated local areas exhibit cooling trends. Although significantly warming areas are discretely distributed, slight cooling trends are present in the central region, yet these do not achieve significance.

4.3 Vegetation change trends based on different vegetation types

Vegetation types in the Qinling–Daba Mountains exhibited pronounced topographic differentiation. Woody vegetation, including broadleaf forests, coniferous forests, and shrubs, dominated mid- to high-altitude mountainous regions, whereas cultivated vegetation and grasslands were concentrated in low-altitude valleys and gentle-slope areas. Cultivated vegetation, broadleaf forests, and shrubs constituted the dominant vegetation types across the region. Forest ecosystems in medium- and high-elevation zones serve as an important ecological barrier for regional ecological stability, while alpine vegetation and marshes are distributed in fragmented patches under hydrologically constrained environments. All vegetation types exhibited positive median Sen’s slopes during 2000–2022 (Figure 9). Cultivated vegetation and grasslands showed the highest growth rates, with median slopes of 0.0039 and 0.0037, respectively, and exhibited widespread high-growth patches. Broadleaf forests showed the lowest median growth rate and a relatively narrow distribution range. Coniferous forests exhibited a pronounced boxplot span and long-tail distribution, indicating substantial differences in vegetation recovery under varying site conditions.

4.4 Trends in vegetation change in the Qinling–Daba Mountains over the past 23 years under different topographic factors

4.4.1 Characteristics of vegetation change at different elevations

Vegetation dynamics in the Qinling–Daba Mountains exhibited pronounced elevational differentiation during 2000–2022 (Figure 10a). NDVI increased across all elevation ranges, with the strongest greening occurring in mid-elevation zones (1,000–2,500 m). The Sen slope peaked at 2,000–2,250 m, indicating the highest vegetation recovery intensity within this range. In contrast, low-altitude (<750 m) and high-altitude (>3,500 m) regions exhibited relatively weak greening trends. Vegetation change exhibited a clear elevation-dependent pattern characterized by degradation at low elevations, rapid greening at mid elevations, and relatively stable dynamics at high elevations (Figure 10b). Low-altitude areas (<750 m) were dominated by decreasing vegetation trends, whereas mid-elevation zones (750–2,250 m) represented the core greening regions with the highest vegetation recovery intensity. At elevations above 2,250 m, vegetation dynamics became increasingly stable, although localized degradation remained detectable in alpine regions above 2,750 m. These patterns indicate substantial elevational differences in vegetation response across the Qinling–Daba Mountains.

4.4.2 Characteristics of vegetation change at different slope gradients

Vegetation changes in the Qinling–Daba Mountains over the past 23 years present a distinct slope-dependent pattern (Figure 11). NDVI generally rises with increasing slope steepness, peaking on steep slopes (30°–35°). In the moderate slope range (15°–20°), NDVI Sen’s slope reaches the highest value (approximately 0.0041) across all slope intervals. Gentle slopes (0°–5°) are dominated by vegetation degradation, with the extremely significant decrease category being the most prominent. The moderate slope range (15°–20°) serves as the primary area of vegetation restoration, with K values for significant and extremely significant increase categories exceeding 1.03. Moderately steep and steep slopes (>15°) consistently show vegetation greening, reflecting substantial recovery of previously disturbed ecosystems. Even on extremely steep slopes (>35°), vegetation still exhibits clear greening trends. This slope-dependent ecological pattern is characterized by degradation on gentle slopes, rapid restoration on moderate slopes, and stable greening on steep slopes.

4.4.3 Characteristics of vegetation change across different aspects

Vegetation distribution and evolution trends in the Qinling–Daba Mountains over the past 23 years show general consistency across different aspects, with only minor local variations. NDVI values increased significantly across all aspects between 2000 and 2022, with small differences among aspect categories (Figure 12a). Mean Sen’s slope values were positive across all aspects (Figure 12b), with northwest-facing slopes having the highest vegetation improvement rate. Aspect also affects the distribution of vegetation change trend categories (Figure 12c). North-, northeast-, and northwest-facing slopes were dominated by vegetation greening, with no significant increase and significant increase categories prevailing. East-, south-, and west-facing slopes featured prominent localized degradation, with south-facing slopes showing the strongest degradation. Localized degradation occurred on south-facing slopes, but vegetation improvement remained dominant across all aspects.

4.5 Evaluation of the potential role of anthropogenic and climatic factors in vegetation change

4.5.1 Bayesian optimization of the XGBoost model

This study utilized an optimized XGBoost model, which exhibited moderate predictive performance across three datasets (Table 3). For the entire study region, the model achieved R² values of 0.59 for the training dataset and 0.43 for the test dataset, with a test RMSE of 0.11 and a test MAE of 0.08. In the region with significant increases, the training and test R² values were 0.48 and 0.39, respectively, accompanied by a test RMSE of 0.10 and a test MAE of 0.08. For the region with significant decreases, the model produced R² values of 0.66 for the training set and 0.33 for the test set, with a test RMSE of 0.16 and a test MAE of 0.13. The lower R² values obtained from the independent test dataset, compared to the training data, suggest a decline in predictive performance when the model was applied to spatially independent samples. This discrepancy highlights the challenges posed by spatial heterogeneity and the complex relationships between NDVI and environmental conditions. Nevertheless, the positive R² values across the three test datasets indicate that the selected environmental variables still provide some predictive information for annual NDVI. Therefore, these optimized models serve as a solid foundation for subsequent SHAP analysis, which will examine the relative contributions and response patterns of the selected environmental drivers.

Table 3. Predictive performance of the optimized XGBoost regression models.

4.5.2 Analysis of drivers of vegetation change

SHAP analysis revealed distinct differences in the relative importance of environmental drivers across the Qinling-Daba Mountains (Figure 13). In the overall study region (Figure 13a), land use and land cover (LULC) emerged as the most influential predictor, accounting for 42.5% of the total mean absolute SHAP contribution. This was followed by digital elevation model (DEM) at 13.2%, population density (Pop) at 8.2%, potential evapotranspiration (PET) at 8.0%, and precipitation (PRE) at 7.0%. In the region experiencing significant increases (Figure 13b), LULC continued to be the dominant predictor, contributing 31.7%. Population density (11.7%), DEM (11.3%), and PRE (11.1%) exhibited relatively similar contributions. Conversely, in the region with significant decreases (Figure 13c), LULC again had the largest contribution at 43.7%, followed by DEM (9.8%), PRE (7.5%), and Pop (7.5%). The relative contributions of the remaining factors were comparatively lower. Overall, land use type consistently demonstrated the highest relative contribution across the three regions, while the importance of other environmental factors varied according to vegetation change status. In the region of significant increases, population density, elevation, and precipitation were relatively prominent. In contrast, in the region of significant decreases, elevation, precipitation, and population density followed LULC in importance. These findings indicate that the relative significance of environmental drivers varied among different vegetation trend zones across the Qinling–Daba Mountains.

SHAP swarm plots further revealed nonlinear and heterogeneous associations between environmental drivers and NDVI across the Qinling–Daba Mountains (Figure 14). In the entire study region (Figure 14a), land use and land cover (LULC) exhibited a broad range of SHAP values, indicating significant variability in its contribution to NDVI predictions. Higher digital elevation model (DEM) values were generally linked to positive SHAP values, while higher population (Pop) values were more frequently associated with negative SHAP values. In the region of significant increase (Figure 14b), LULC, DEM, and Pop displayed relatively pronounced response patterns. Here, higher DEM values were typically connected to positive contributions, whereas higher Pop values were more often associated with negative contributions. Additionally, precipitation (Pre) and potential evapotranspiration (PET) demonstrated nonlinear variations in their SHAP values. In the region of significant decrease (Figure 14c), LULC again showed considerable variation in SHAP values. Meanwhile, DEM was generally linked to positive contributions, and Pop was associated with negative contributions at higher feature values. Other climatic variables, including Pre, soil moisture (SM), PET, and vapor pressure deficit (VPD), exhibited heterogeneous responses across their respective feature ranges. Overall, these findings indicate that the relationships between environmental conditions and NDVI are nonlinear and vary among different vegetation trend zones.

4.5.3 Univariate analysis of vegetation changes

The land use and land cover (LULC) factor, which has the highest contribution, exerts directional effects on the normalized difference vegetation index (NDVI) that vary by land use type. In the entire study region (Figure 15a), forest consistently demonstrated a positive SHAP contribution to NDVI. In contrast, cropland, shrubland, grassland, and impervious land were primarily associated with negative contributions. Among these land-use types, impervious land exhibited the most pronounced negative contribution, while forest showed the strongest positive contribution. In the region of significant increase (Figure 15b), forest again displayed a consistently positive SHAP contribution. Conversely, cropland, shrubland, and grassland predominantly exhibited negative contributions. Compared to the entire study region, the magnitude of negative contributions varied among these land-use types, yet the positive contribution of forest remained evident. Notably, only four LULC categories were represented in this subset, as impervious land was absent from the test samples in the significant increase region. In the region of significant decrease (Figure 15c), forest again showed the strongest positive SHAP contribution. Cropland and shrubland also exhibited predominantly positive contributions. In contrast, grassland displayed a predominantly negative contribution, while impervious land had the strongest negative contribution among the five land-use categories. These results indicate that SHAP contributions associated with individual land-use types varied across vegetation trend zones, particularly for cropland, shrubland, and grassland.

The SHAP dependence plots demonstrated nonlinear relationships between the driving factors and NDVI, revealing distinct response patterns across the overall, significant increase, and significant decrease regions (Figure 16). For the Digital Elevation Model (DEM), SHAP contributions generally shifted from negative to positive as elevation increased in the overall region (Figure 16a). In the significant increase region, contributions rose sharply at lower elevations, reached positive values around 900 m, and remained relatively stable thereafter. Conversely, the significant decrease region displayed a more intricate non-monotonic pattern. Here, contributions initially declined from positive to negative at intermediate elevations before increasing again at higher elevations. PRE exhibited a shift from negative to positive SHAP contributions as precipitation increased in both the overall and significant increase regions (Figure 16b). In contrast, the significant decrease region displayed a non-monotonic response. Here, contributions were relatively high at low precipitation levels, decreased to near-zero or negative values at intermediate levels, and gradually increased at higher precipitation. Pop primarily showed decreasing SHAP contributions with increasing density in both the overall and significant increase regions (Figure 16c). Contributions turned negative at approximately 2.1–2.2 (~14 people/km²). In the significant decrease region, contributions remained positive over a broader density range before declining at higher densities. SM demonstrated contrasting response patterns (Figure 16d). In the significant increase region, rising SM reduced SHAP contributions from positive to zero near the 0.35 threshold, becoming negative above this value. Conversely, in the significant decrease region, SHAP contributions increased from negative to positive once SM exceeded 0.31 and stabilized at higher SM values. Overall, SM exerted a weak effect throughout the study; SHAP contributions remained near zero across most SM ranges, with the transition from positive to negative effects occurring around 0.38. These findings indicate that the modeled relationships between environmental and anthropogenic factors and NDVI were nonlinear and varied significantly across different vegetation trend regimes.

5 Discussion

5.1 Elevational differentiation of vegetation recovery

Vegetation recovery in the Qinling–Daba Mountains exhibits variation with altitude, peaking at mid-levels. At lower elevations, human activities impede recovery, while higher altitudes are limited by thermal constraints. This pattern results from the interplay between hydrothermal conditions and anthropogenic disruptions across the landscape. Although lower elevations possess favorable hydrothermal conditions, they are subject to significant human pressures. Consequently, vegetation recovery is restricted, leading to localized degradation (Xu J. et al., 2020). In contrast, mid-levels serve as transition zones with optimal ecological conditions. Here, reduced human disturbances and effective restoration efforts promote vegetation growth (Cheng et al., 2021; Liu H. et al., 2024). Meanwhile, higher altitudes experience minimal human influence; however, they face challenges such as cold temperatures and limited thermal resources, which hinder vegetation growth (Liu et al., 2018).

Despite overall greening trends, recovery rates and stability are lower at high elevations compared to mid-elevation zones. Variations in elevational adaptability among different vegetation types further exacerbate the spatial heterogeneity of recovery. Shallow-rooted vegetation, including grasses and crops, is more tolerant of low- and mid-elevation conditions and exhibits faster recovery (Tuo et al., 2024). In contrast, forest ecosystems, characterized by their complex structure and longer regeneration cycles, demonstrate more spatially heterogeneous responses. In summary, vegetation recovery in mountainous regions is influenced by hydrothermal suitability and the elevational distribution of human disturbance.

5.2 Climate–water constraints on vegetation dynamics

Climatic factors, especially atmospheric water stress, significantly influence vegetation dynamics in the Qinling–Daba Mountains. Nonlinear threshold effects characterize vegetation responses to climate, with ecological thresholds playing a critical role in regional vegetation degradation and fluctuation. The prevailing warm–humid climate trend in this area fosters vegetation greening by extending the growing season and enhancing photosynthetic efficiency (Yao and Cui, 2022). Sufficient precipitation replenishes soil moisture, thereby promoting growth. However, surpassing critical thresholds can yield detrimental effects. For instance, while an improved water supply can mitigate drought stress and enhance photosynthesis and leaf growth when soil moisture is low, exceeding a critical threshold can impair soil aeration. This impairment may result in anaerobic waterlogging, which damages the root system and inhibits vegetation growth (Kong et al., 2017). Furthermore, increased evapotranspiration and vapor pressure deficit can create “warming without wetting” conditions (Guo et al., 2021; Yuan et al., 2019). The benefits of precipitation are limited by thermal constraints, which stabilize vegetation growth. This response, shaped by hydrothermal imbalance and atmospheric water stress, elucidates localized vegetation degradation in warm-humid climates.

5.3 Human activities and land-use regulation

Human activities and land-use changes significantly influence vegetation dynamics in the Qinling–Daba Mountains. Vegetation recovery demonstrates nonlinear responses to varying levels of human disturbance. Intensive land-use practices, such as urban expansion and agriculture, disrupt natural land covers and lead to landscape fragmentation (Liu H. et al., 2024). This fragmentation impedes vegetation recovery, particularly in low-elevation valleys and urban peripheries.

Ecological restoration initiatives, including the Grain-for-Green Program and forest conservation efforts, have successfully promoted vegetation recovery in mid-elevation mountainous regions (Yao, 2024). Population density affects vegetation dynamics in two ways. In high-density areas, increased development pressure results in vegetation degradation (Chen et al., 2025). In contrast, regions with lower human disturbance and moderate management practices facilitate vegetation recovery. The ecological effects of human activities vary according to disturbance levels and land use. Vegetation changes in the Qinling–Daba Mountains are shaped by environmental factors, shifts in land use, and human activities (Yao and Cui, 2022). Therefore, controlling the intensity of development and maintaining ecological restoration efforts are essential for the stability of mountain ecosystems.

5.4 Research uncertainties and future perspectives

This study examined vegetation dynamics and their climatic, topographic, and anthropogenic drivers in the Qinling–Daba Mountains utilizing long-term NDVI data and a machine-learning framework. Several uncertainties warrant acknowledgment. First, NDVI primarily indicates vegetation greenness and may saturate in densely vegetated areas. This limitation restricts its ability to represent biomass and ecosystem functioning (Kang et al., 2021). Additionally, annual maximum NDVI fails to capture seasonal dynamics or time-lagged vegetation responses. Second, aggregating multi-source datasets to a 1-km analysis scale may obscure fine-scale heterogeneity, particularly in 30-m NDVI and topographic variables, thereby introducing scale-dependent uncertainty. Third, multicollinearity among predictors, especially DEM, TMP, PET, SM, and VPD, may influence SHAP-based importance rankings. Consequently, the identified relationships should be viewed as model-based associations rather than independent causal effects. Fourth, although a 10-km spatial-block validation was implemented to mitigate spatial dependence and information leakage, the original random sampling design may still restrict sample representativeness across environmental gradients. Furthermore, factors such as ecological restoration and land management were not explicitly considered (Shi et al., 2020).

Future studies should integrate multi-seasonal vegetation indices and lagged environmental variables to enhance the characterization of delayed vegetation responses. Additionally, higher-resolution observations and field measurements may help mitigate scale-related uncertainties. Employing spatially stratified sampling along with more rigorous spatial cross-validation could enhance sample representativeness and model robustness. Furthermore, alternative attribution methods or grouped predictors could be utilized to evaluate the stability of SHAP-based importance rankings. Including detailed information on ecological restoration and human activities would also deepen the understanding of vegetation dynamics in the context of environmental change.

6 Conclusion

This study examined vegetation dynamics and their potential climatic, topographic, and anthropogenic drivers in the Qinling–Daba Mountains, utilizing 23 years of annual maximum NDVI data from 2000 to 2022. The primary conclusions are as follows:

• The Qinling–Daba Mountains experienced a significant overall greening trend during the period from 2000 to 2022, with NDVI increasing at approximately 0.0051 yr.⁻¹ (p < 0.05). Approximately 90.13% of the study area exhibited an increasing NDVI trend, and Hurst analysis indicated that 89.83% of the area is likely to sustain this improving trend. Notably, vegetation improvement was most pronounced in regions with moderate elevations (1,750–2,250 m), moderate slopes (15–20°), and north-facing slopes. In contrast, stronger human influences were primarily observed in low-elevation and steep-slope areas.
• Topographic heterogeneity contributed to significant spatial variation in vegetation conditions. Elevation demonstrated a more substantial relationship with NDVI compared to slope and aspect, while the influences of slope and aspect were relatively weaker. These findings suggest that vegetation dynamics varied along topographic gradients, underscoring the necessity of considering terrain heterogeneity when assessing vegetation changes in mountainous ecosystems.
• SHAP analysis identified land-use type (LULC) as the most influential predictor across all three regions, contributing 42.5, 31.7, and 43.7% of the mean absolute SHAP values in the overall, significant increase, and significant decrease regions, respectively. In the overall study area, the next most important predictors were elevation (13.2%), population density (8.2%), potential evapotranspiration (PET) (8.0%), and precipitation (7.0%). The relative importance of these factors varied across vegetation trend zones, indicating heterogeneous associations between NDVI and environmental conditions. SHAP dependence analysis further revealed nonlinear response patterns, including an approximately 850-m transition in the elevation–NDVI relationship in both the overall and significant increase regions. Additionally, there were nonlinear responses to precipitation, population density, and soil moisture.
• Overall, vegetation dynamics in the Qinling–Daba Mountains were characterized by widespread greening, yet responses varied across environmental gradients. These results emphasize the importance of jointly considering land-use conditions, elevation, climatic variability, and human activity when assessing vegetation dynamics in ecological transition zones. Such findings can inform spatially differentiated ecological restoration and vegetation management. However, the identified model-based associations should be interpreted cautiously, rather than as direct causal relationships.

Generative AI statement

The author(s) declared that Generative AI was used in the creation of this manuscript. During the preparation of this manuscript, the authors used DeepL Write and ChatGPT 4o to assist with language polishing and grammatical correction. The tool was used solely to improve the clarity and readability of the text, without any modification to the core research content, data interpretation, or scientific conclusions. After using the tool, the authors reviewed and edited the entire content and take full responsibility for the accuracy and integrity of the published work.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Petras Rupšys, Vytautas Magnus University, Lithuania

Reviewed by: Adil Dilawar, Beijing Normal University, China

Zihe Wang, China Institute of Water Resources and Hydropower Research, China