Trend Analysis and Driving Factors of Vegetation Dynamics in Northern China from 1982 to 2015

Abstract

Under the background of global warming, understanding the dynamic of vegetation plays a key role in revealing the structure and function of an ecosystem. Assessing the impact of climate change and human activities on vegetation dynamics is crucial for policy formulation and ecological protection. Based on the Global Inventory Monitoring and Modeling System (GIMMS) third generation of Normalized Difference Vegetation Index (NDVI3g), meteorological data and land cover data, this study analyzed the linear and nonlinear trends of vegetation in northern China from 1982 to 2015, and quantified the relative impact of climate change and human activities on vegetation change. The results showed that more than 53% of the vegetation had changed significantly, and 36.64% of the vegetation had a reverse trend. There were potential risks of vegetation degradation in the southwestern, northwestern and northeastern parts of the study’s area. The linear analysis method cannot disclose the reversal of the vegetation growth trend, which will underestimate or overestimate the risk of vegetation degradation or restoration. Climate change and human activities promoted 76.54% of the vegetation growth in the study area, with an average contribution rate of 51.22% and 48.78%, respectively, while the average contribution rate to the vegetation degradation area was 47.43% and 52.57%, respectively. Vegetation restoration of grassland and woodland was mainly affected by climate change, and human activities dominated their degradation, while cropland vegetation was opposite. The contribution rate of human activities to vegetation change in the southeastern and eastern parts of the study area was generally higher than that of climate change, but it was the opposite in the high altitude area, with obvious spatial heterogeneity. These results are helpful to understand the dynamic mechanism of vegetation in northern China, and provide a scientific basis for vegetation restoration and protection of regional ecosystems.

1. Introduction

The dynamic change of vegetation can affect global carbon, the water cycle and energy flow [1], and is considered an indicator of ecosystems’ health [2,3]. Vegetation plays an indispensable role in maintaining regional ecological balance, especially in reducing soil erosion [4]. Therefore, monitoring vegetation dynamics and analyzing their drivers are very important for studying global climate change and ecological environment protection and restoration [5,6].

The rapid development of satellite remote sensing provides a series of continuous spatiotemporal observation data for scientific research [7]. Remote sensing data have become an important source in the field of ecological environment monitoring because of their advantages such as large space coverage, easy access, rich information and continuous time scale [8,9]. The normalized difference vegetation index (NDVI) calculated from red and near-infrared reflectances is considered an effective indicator to evaluate vegetation growth state [7], and has good correlation with leaf area, photosynthetically active radiation absorption, and vegetation productivity [10]; therefore, it is often used to study vegetation dynamics and monitor the ecological environment [5,11]. The increase (greening) or decrease (browning) of NDVI can directly reflect the enhancing or weakening of vegetation activity [12]. Therefore, correctly describing vegetation dynamics is the key to understanding the impact of climate change on terrestrial ecosystems [12].

Northern China includes several typical areas, such as the Loess Plateau, Xinjiang, Inner Mongolia, Qilian Mountains and Qinghai Plateau, which are the key areas for water and soil conservation, vegetation protection and restoration [13]. In the past few decades, in order to protect and improve the regional ecological environment, the Chinese government has launched a series of large-scale ecological protection and restoration projects, such as the “Grain for Green Program” [13], “Natural Forest Conservation Program” [14], “Three North Shelterbelt Development Program” [15], and the “Beijing-Tianjin Sand Source Control Program” [16]. The implementation of ecological projects can significantly improve the vegetation coverage and ecological environment [17,18]. The driving factors of vegetation dynamics can be divided into climate factors and human activities. An accurate description of the vegetation dynamic process is helpful to accurately evaluate the relative contribution of climate change and human activities to vegetation dynamics, and to divide the different driving forces of vegetation dynamics, which is of great significance for ecosystem protection and policy formulation [7].

At present, simple linear models are used by many researchers to describe the process of vegetation change, assuming that the process of vegetation change is monotonic with a constant trend [12]. However, the growth process of vegetation is unstable and nonlinear, which changes with time [19]. Therefore, the linear trend analysis method may ignore the potential greening or browning trend of vegetation during the study period [10,12]. A nonlinear analysis method is considered effective for analyzing vegetation dynamics [7], mainly including piecewise linear regression model [20], breaks for additive seasonal and trend (BFAST) [21], detecting breakpoints and estimating segments in trend (DBEST) [22], etc. The above methods are very sensitive to short-term fluctuations and abrupt changes, which weakens their ability to detect trend shifts and turning points in long time series [12]. Wavelet analysis needs to select basic functions, which also limits its application in vegetation dynamics [23]. The ensemble empirical mode decomposition (EEMD) [24] is usually used to study the time series trend of variables in signal and image processing [25], climate diagnosis [26], hydrology [27] and other fields. EEMD does not need to determine the basic function, and has good adaptability [12]. At the same time, EEMD decomposition can reveal more potential information of nonlinear and non-stationary time series of vegetation [6], which has been applied in the remote sensing field [10,23,28,29].

The vegetation types in northern China are complex and diverse, with most of them being located in semi-arid and arid areas, and the natural ecological environment is very complex and fragile [30]. Therefore, analyzing the nonlinear change characteristics of vegetation will help us to better understand the dynamics of vegetation and its relationship with driving factors [7]. At present, most of the relevant studies focus on typical ecological study areas [7,31], while few studies focus on the whole of northern China, lacking a comprehensive analysis of the spatial distribution of vegetation greening, browning and trend reversal. Therefore, in order to deeply understand the vegetation dynamics and its driving mechanism in northern China, the purpose of this study is to (a) analyze the spatial pattern of linear and nonlinear changes in vegetation, (b) analyze the nonlinear characteristics of vegetation in different land cover and climatic zones, (c) assess the response of vegetation with different trends to climate change, and (d) quantify the relative contributions of climate change and human activities to vegetation change.

2. Materials and Methods

2.1. Study Area

The study area is located in northern China, which lies between 73°37′–125°57′E and 31°41′–53°18′N, with a total area of about 4,497,784 km². This study area covers three steps of China’s terrain, with obvious differences in altitude, ranging from −156 to 7347 m. The altitude in the west and southwest of the study area is relatively high, while that in the east is relatively low. The monthly average temperature is about 4.37 °C, and the monthly average precipitation is about 23.48 mm. According to the precipitation, it can be divided into arid, sub-arid, sub-humid and humid areas, of which arid and semi-arid areas are the main types of the study area. According to the climatic zone, it is divided into subtropical, warm temperate, middle temperate, cold temperate, and plateau climate zones. The study area is mainly covered by warm temperate, middle temperate and plateau climate zones [32] (Figure 1). Therefore, there are regional differences in hydrothermal environment in the study area (

Table 1. Detailed information on 11 climate zones in the study area.

Code	Climate Condition	Code	Climate Condition
CTH	Cold-temperate humid zone	WTA	Warm temperate arid zone
MTH	Middle temperate humid zone	PCSH	Plateau climate sub-humid zone
MTSH	Middle temperate sub-humid zone	PCSA	Plateau climate sub-arid zone
MTSA	Middle temperate sub-arid zone	PCA	Plateau climate arid zone
MTA	Middle temperate arid zone	SH	Subtropical humid zone
WTSH	Warm temperate sub-humid zone

), forming a unique terrestrial ecosystem. Grassland is the main land cover type in the study area (Figure 1c).

2.2. Data and Preprocessing

2.2.1. GIMMS NDVI3g

GIMMS NDVI3g, with 8 km spatial resolution and 15-day temporal resolution, is composed of upper and lower scene data every month, and the time span ranges from July 1981 to December 2015. Although the spatial resolution of GIMMS NDVI3g is relatively low, its time series is the longest and can reflect the change trend of NDVI more comprehensively than other products [33]. Therefore, GIMMS NDVI3g has been widely used for global or regional scale vegetation change monitoring [12]. High quality NDVI time series data are of great significance for regional and global ecological and environmental application research. In this study, spatial–temporal Savitzky–Golay (STSG) filter proposed by Cao [34] was used to smooth NDVI data, so as to further improve the quality of GIMMS NDVI3g. Finally, the filtered GIMMS NDVI3g was aggregated into monthly NDVI data by Maximum Value Composite (MVC) [35]. Combined with previous studies, the vegetation growth season was defined as April to October [36], and the NDVI in the growth season (GSN) was obtained by calculating the average value from April to October.

2.2.2. Meteorological Datasets

The meteorological data included precipitation and temperature, which came from the daily value dataset of surface climate data (V3.0) provided by China Meteorological Data Service Center (http://data.cma.cn/ (accessed on 21 May 2018)), including 829 meteorological stations, with a time range from 1981 to 2015. The monthly raster data with the same spatial resolution as GIMMS NDVI3g were obtained by interpolating the meteorological station data with digital elevation model (DEM) data and ANUSPLIN 4.2 software [37]. The DEM data were collected from the NASA Shuttle Radar Topographic Mission [23].

The monthly radiation raster data were downloaded from ERA5 of the European Centre for Medium-Range Weather Forecasts (ECMWF) (https://www.ecmwf.int/ (accessed on 12 July 2019)), which was the fifth generation ECMWF atmospheric reanalysis of the global climate covering the period from January 1950 to present, with a spatial resolution of 0.25 degrees. Similarly, we used ANUSPLIN 4.2 software and DEM data to interpolate the solar radiation data to match the spatial resolution of GIMMS NDVI3g.

Finally, we calculated the mean temperature (GS-Tem), cumulative precipitation (GS-Pre) and mean solar radiation (GS-SR) in the growing season of vegetation.

2.2.3. Other Geospatial Ancillary Data

The land cover data and climate zone data were derived from Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC, http://www.resdc.cn (accessed on 1 August 2019)). The land cover data contain seven periods of land cover data, with a spatial resolution of 1 km and a time span from 1980 to 2015. The first-class classification accuracy exceeds 94.3% [38], and has been widely used in relevant research on a regional scale [30,39,40]. The land cover data were resampled to the same spatial resolution as GIMMS NDVI3g using the nearest neighbor resampling method. In order to reduce the impact of land cover change and classification errors, we only selected the grids with unchanged land cover types from 1980 to 2015. The vegetation cover types included cropland, woodland and grassland, which accounted for 18.43%, 16.02% and 65.55% of the total area (referring to the total area of cropland, woodland and grassland).

2.3. Method

2.3.1. Sen Median Trend Analysis and Mann–Kendall Test

Sen median trend analysis and the Mann–Kendall test are widely used in trend analysis of meteorological, hydrological and remote sensing time series [33,41]. Sen median trend analysis can avoid the interference of outliers well, which is more suitable for studying vegetation change trend [42]. The calculation formula can refer to previous studies [42].

The Mann–Kendall test was used to judge the significance of time series trend. The formula of the Mann–Kendall test can refer to previous studies [30]. In the study, $\alpha$ was set to 0.05, and the NDVI trend was divided into three types: monotonic browning (B to B) (ρ < 0, |Z| > 1.96, p < 0.05), monotonic greening (G to G) (ρ > 0, |Z| > 1.96, p < 0.05) and non-significant change (Non-sig) (|Z| ≤ 1.96, p > 0.05).

2.3.2. Ensemble Empirical Mode Decomposition (EEMD) Method

Wu [24], aiming at overcoming the shortcomings of empirical mode decomposition (EMD), proposed ensemble empirical mode decomposition (EEMD) based on EMD. EEMD uses Gaussian white noise to assist the original data. After the amplitude of white noise meets certain conditions, after several calculations, the set average is enough to make the white noise cancel each other, and its value will not have a significant impact on the decomposition results. EEMD method has strong adaptability and avoids human experience interference [12]. Therefore, we used EEMD to study the nonlinear change characteristics of vegetation, and obtained the trend component that could reflect the long-term trend of the original time series data. The decomposition steps of EEDM can refer to previous studies [12,29].

Considering the calculation time and robustness of EEMD [12], the number of Gaussian white noise added in this study was set to 100, and the noise amplitude was 0.2 standard deviation of the original time series data.

In order to compare with the linear trend, we used the average changing rate of instantaneous trend as the average EEMD trend of NDVI time series during the study period [12,28], which was the EEMD trend of the last year divided by the time interval.

$${{\displaystyle T\mathrm{rend}}}_{eemd} = \frac{{{\displaystyle R}}_n({{\displaystyle t}}_{end})-{{\displaystyle R}}_n({{\displaystyle t}}_0)}{{{\displaystyle t}}_{end}-{{\displaystyle t}}_0}$$

(1)

where t₀ and t_end represent the time at the beginning and end of the study, respectively.

Similarly, we used the method of Pan et al. [12] to test the significance of EEMD trend. According to the significance test and extreme points, the EEMD trend of NDVI was divided into five types:

(1)

Non-significant (Non-sig): the trends were not significant at any year (p > 0.05).

(2)

Greening to greening (G to G): the trends were monotonic increasing and were statistically significant for at least one year (p < 0.05).

(3)

Browning to browning (B to B): the trends were monotonic decreasing and were statistically significant for at least one year (p < 0.05).

(4)

Greening to browning (G to B): the trends contained one local maximum, which changed from increasing trend to decreasing trend, and were statistically significant for at least one year (p < 0.05).

(5)

Browning to greening (B to G): the trends contained one local minimum, which changed from decreasing trend to increasing trend, and were statistically significant for at least one year (p < 0.05).

Four pixels in the study area were randomly selected to show different types of NDVI EEMD trends (Figure 2). The trend component R_n represents the long-term change trend of NDVI time series. For trends with one extreme value, the location of the extreme value is identified as the turning point (TP). Monotonic greening and browning were collectively referred to as monotonic trend. Greening to browning and browning to greening were collectively referred to as reversal trend.

2.3.3. Partial Correlation Analysis

For multivariate correlation analysis, Pearson correlation analysis cannot exclude the impact of mutual correlation between other climate factors on NDVI. Therefore, a partial correlation coefficient (PCC) was used to calculate the impact of a single climate factor on vegetation growth by treating other climate factors as control variables. For three climate factors, PCC is as follows:

$${{\displaystyle R}}_{12,34} = \frac{{{\displaystyle R}}_{12,3}-{{\displaystyle R}}_{14,3}\times {{\displaystyle R}}_{24,3}}{\sqrt{(1-{{\displaystyle R}}_{14,3}^2)\times (1-{{\displaystyle R}}_{24,3}^2)}}$$

(2)

$${{\displaystyle R}}_{12,3} = \frac{{{\displaystyle R}}_{12}-{{\displaystyle R}}_{13}\times {{\displaystyle R}}_{23}}{\sqrt{(1-{{\displaystyle R}}_{13}^2)\times (1-{{\displaystyle R}}_{23}^2)}}$$

(3)

where $R_{12,34}$ refers to the PCC of variables 1 and 2, and variables 3 and 4 are the control variables;${ R}_{12,3}$ refers to the PCC of variables 1 and 2, and variable 3 is the control variable;${ R}_{12}$ refers to the correlation coefficient of variables 1 and 2. Similarly, other variables have similar meanings. Variable 1 refers to NDVI, and variables 2, 3, and 4 represent three climatic factors, respectively. The PCC between NDVI and each climatic factor was calculated. The t-test was used to estimate the significance of the calculated PCC in the case of p < 0.05.

2.3.4. Residual Trend (RESTREND) Analysis

The driving factors of vegetation dynamics include climate and non-climate factors, with non-climate factors usually being considered as human activities [43]. The residual trend analysis method is usually used to distinguish the impact of human activities and climate change on vegetation change [39,44]. Firstly, multiple linear regression is used to simulate the predicted value of NDVI under the influence of climate change (NDVI_c). Then, calculate the NDVI value affected by human activities (NDVI_h), which is the difference between the actual NDVI and NDVI_c. The positive or negative slopes of NDVI_c and NDVI_h indicate whether climate change and human activities promote vegetation growth, respectively.

$${\mathrm{NDVI}}_c = p_0 + p_{1}Tem + p_{2}Pre + p_{3}SR + \epsilon$$

(4)

$${\mathrm{NDVI}}_h = \mathrm{NDVI} - {\mathrm{NDVI}}_c$$

(5)

where Tem, Pre, and SR are the temperature, precipitation, and solar radiation, respectively; $p_0$ $p_1$, $p_2$ and $p_3$ are the fitting coefficients; and ε is the residual error term.

According to previous studies [28,41,44], climate change and human activities can be divided into six scenarios according to the positive and negative relationship among NDVI, NDVI_c and NDVI_h slopes, as shown in

Table 2. Assess the relative contribution of climate change (CC) and human activity (HA) drivers in vegetation dynamics under six scenarios.

SNDVI	SCC	SHA	Relative Contribution of CC	Relative Contribution of HA	Driving Forces of Vegetation Dynamics
+	+	+	$\frac{\left|S_{CC}\right|}{\left|S_{CC}\right|+\left|S_{HA}\right|}$	$\frac{\left|S_{HA}\right|}{\left|S_{CC}\right|+\left|S_{HA}\right|}$	Both CC and HA induced NDVI increase (ICH)
	+	−	100	0	CC induced NDVI increase (ICC)
	−	+	0	100	HA induced NDVI increase (IHA)
−	−	−	$\frac{\left|S_{CC}\right|}{\left|S_{CC}\right|+\left|S_{HA}\right|}$	$\frac{\left|S_{HA}\right|}{\left|S_{CC}\right|+\left|S_{HA}\right|}$	Both CC and HA induced NDVI decrease (DCH)
	−	+	100	0	CC and induced NDVI decrease (DCC)
	+	−	0	100	HA induced NDVI decrease (DHA)

Abbreviations: S_NDVI, S_CC and S_HA represent the slopes of NDVI, NDVI_c and NDVI_h, respectively.

.

To better explain the overall process of this study, Figure 3 shows the flow chart.

3. Result

3.1. Temporal Variations of GSN from 1982 to 2015

The linear trend showed that the average GSN fluctuated and increased significantly, with an average rate of 0.008/10a (p < 0.01) (Figure 4). Among them, the increasing rate of cropland GSN was the highest and that of grassland vegetation was the lowest, which were 0.016/10a (p < 0.01) and 0.006/10a (p < 0.01), respectively. According to EEMD method, GSN of all vegetation land also showed a significant increase trend (p < 0.05), and the change rate increased over time from 0.010/10a in 1983 to 0.013/10a in 2015. The EEMD trend of GSN of the three vegetation types increased significantly (p < 0.05), and the change rates of cropland and grassland increased over time, from 0.018/10a and 0.007/10a in 1983 to 0.020/10a and 0.012/10a in 2015, respectively, while the change rate of woodland decreased from 0.007/10a in 1983 to 0.006/10a in 2015. In order to compare with the linear trend, we calculated their average EEMD trend, which were 0.012/10a, 0.019/10a, 0.006/10a and 0.009/10a, respectively. The average EEMD trend of cropland and grassland was significantly higher than the linear trend, while that of woodland was the opposite. The EEMD trends of cropland and grassland in the later stage were higher than those in the earlier stage, indicating vegetation greenness accelerating, while the EEMD trend of woodland was decreasing, indicating slow vegetation greenness.

3.2. Spatial Pattern of Linear and Nonlinear Trends of GSN

Figure 5a,b show that the spatial distribution of high and low values of the change rate of Sen trend and average EEMD trend is relatively similar. The linear trend indicated that the area with a significant monotonic greening trend (G to G) was slightly larger than that with a non-significant trend (Non-sig), accounting for 47.69% and 46.55% of the total area, respectively. The area with a significant monotonic browning trend (B to B) was the lowest, accounting for only 5.76% (Figure 5c). The regions with a monotonic greening trend of GSN were mainly distributed in the WTSH, SH, PCA, and the southern MTA. The areas with a monotonic browning trend of vegetation were mainly concentrated in the northeastern and southwestern parts of the study area.

The area with vegetation monotonic trend (monotonic greening and monotonic browning) detected by EEMD accounted for 23.41% of the total area, of which the area with monotonic greening trend accounted for 22.29%, and the area with monotonic browning trend accounted for only 1.12%. According to EEMD, 36.64% of the region’s vegetation growth trend had shifted, of which 16.17% had changed from greening to browning, and 20.47% had changed from browning to greening. The areas with vegetation monotonic greening trend were mainly distributed in the southeastern and southern parts of the study area, while monotonic browning areas were small and scattered. The areas with greening to browning trend of vegetation mainly occurred in the northeastern, eastern, northwestern and southwestern parts of the study area. The areas with browning to greening trend were mainly distributed in the central and western parts of the study area and the northern of the plateau climate zone (Figure 5d).

Regardless of linear or nonlinear methods, the area with vegetation monotonic browning trend was the smallest, accounting for 1.12% and 5.76% of the total area, respectively (

Table 3. Superposition statistics of linear and nonlinear trend type results (%).

Trend Types		Nonlinear Trend Base on EEMD
		Non-Sig	G to B	B to G	G to G	B to B	Total
Mann-Kendall result	B to B	2.02	2.65	0.19	0	0.90	5.76
	Non-sig	30.70	9.36	4.50	1.77	0.22	46.55
	G to G	7.23	4.16	15.78	20.52	0	47.69
	Total	39.95	16.17	20.47	22.29	1.12	100

Abbreviations: Non-sig: Non-significant, G to B: Greening to browning, B to G: Browning to greening, G to G: Monotonic greening, B to B: Monotonic browning.

). Compared with EEMD method, the Mann–Kendall method can only distinguish three types (B to B, G to G and Non-sig), and it is impossible to judge whether there is a reversal trend of vegetation growth. For example, in the central and eastern regions of the study area, the results of the Mann–Kendall method showed a monotonic greening, while the trend obtained by EEMD method was from browning to greening and from greening to browning. The vegetation in about 30.70% of the study area showed a non-significant trend using the two methods, but the vegetation in 9.36%, 4.50%, 1.77% and 0.22% of the areas showed greening to browning trend, browning to greening trend, monotonic greening trend and monotonic browning trend for the EEMD method, while the Mann–Kendall method showed a non-significant change. About 20.52% of the vegetation in the study area showed a monotonic green trend for the two methods, but about 7.23%, 4.16% and 15.78% of the vegetation showed a non-significant trend, greening to browning trend and browning to greening trend for the EEMD method, while the Mann–Kendall method showed a significant monotonic greening trend. Only 0.90% of the vegetation in the study area showed monotonic browning trend for the two methods, but about 2.02%, 2.65% and 0.19% of the vegetation for the EEMD method was a non-significant trend, greening to browning trend and browning to greening trend, while the Mann–Kendall method showed a significant monotonic browning trend. The area percentage of monotonic greening trend for the Mann–Kendall method was much larger than that for the EEMD method, which was 47.69% and 22.29%, respectively (Table 3). The results showed that EEMD is necessary to analyze the actual evolution of vegetation and the deficiency of linear trend, while linear trend analysis underestimates the potential possibility of vegetation restoration or degradation.

Based on EEMD, the TP of vegetation growth was obtained, which mainly occurred before 1995, accounting for 63.65% of the area with the reversal trend, especially before 1990, which accounted for 38.29%, mainly concentrated in the central part of the study area (Figure S1). The average EEMD trend of vegetation increased from 0.004/10a before the TP to 0.008/10a after the TP. The trend increased significantly, which was twice more than before the TP. Before the TP, the trend of vegetation degradation was mostly concentrated in the range of −0.001 to 0 (Figure 6). Before and after the TP, the browning to greening trend can be divided into monotonic browning and monotonic greening trend, and the greening to browning trend can be divided into monotonic greening and monotonic browning trend. The vegetation restoration trend was obvious, from −0.007/10a of browning trend to 0.026/10a of greening trend. Similarly, the risk of vegetation degradation was very high, from 0.018/10a of greening trend to −0.016a/10a of browning trend. After the TP, the vegetation growth and restoration effect in the central part of the study area was obvious and maintains a high greening trend, while the vegetation changed from a high greening trend to a high browning trend in the southwestern, northwestern, and northeastern parts of the study area, with a high risk of degradation (Figure 6). Therefore, it is necessary to pay more attention to the shift trend of vegetation and evaluate the degradation risk of different vegetation.

3.3. EEMD Trend of Vegetation in Different Climatic Zones and Land Cover

The cropland area with an increasing trend (monotonic greening and browning to greening) accounted for 56.67%, which was much higher than that with degradation or potential degradation trend (monotonic browning or greening to browning) (accounting for 19.20%) (

Table 4. The trend of vegetation GSN in different land cover types for EEMD (%).

Land Cover Types	Nonlinear Trend Base on EEMD
	Non-Sig	G to B	B to G	G to G	B to B	Total
Cropland	24.13	18.28	22.77	33.90	0.92	100
Woodland	42.44	20.63	12.61	23.06	1.26	100
Grassland	43.80	14.48	21.74	18.84	1.14	100

Abbreviations: Non-sig: Non-significant, G to B: Greening to browning, B to G: Browning to greening, G to G: Monotonic greening, B to B: Monotonic browning.

). The vegetation of many croplands had changed from previous degradation to restoration. The trend of woodland was mainly non-significant, accounting for 42.44% in area. Among the significant trends, the monotonic greening area was the largest (23.06%), followed by greening to browning (20.63%). Similarly, the area of grassland with the non-significant trend accounted for 43.80%. Among the significant trends, the browning to greening area accounted for the largest (21.74%), followed by monotonic greening (18.84%). The areas of cropland, woodland and grassland with reverse trend were higher than that with monotonic trend, which accounted for 41.05% vs. 34.82%, 33.24% vs. 24.32%, 36.22% vs. 19.98%, respectively. The potential degradation risk of woodland was higher, followed by cropland (Table 4).

The cropland in SH, WTSH, WTA and MTSH showed an obvious increasing trend, accounting for 73.06%, 61.58%, 62.33% and 57.25% of the total cropland area in the corresponding climate zones, respectively (Table S1). The increasing trend of woodland in SH and WTSH was obvious, accounting for 54.70% and 63.54% of the total woodland area in the corresponding climate zone, respectively, of which the monotonic greening trend accounted for 44.16% and 44.52%, respectively. The area proportion of the woodland with the non-significant trend in CTH, MTSH and MTSA was higher than 50%, and the potential risk of vegetation degradation was high. The potential degradation risk of grassland in PCSH was the highest, followed by the CTH (the number of pixels in MTH was small and could be ignored), and the areas from greening to browning trend were 48.49% and 30.68%, respectively. Grassland had been significantly improved in MTA, WTSH, WTA, PCA and SH, accounting for 43.60%, 79.92%, 42.13%, 58.83% and 73.94% of the grassland area in the corresponding climate zones, respectively. Among them, the monotonic greening trend of grassland in WTSH and SH was 42.06% and 58.79, respectively, which was significantly better than that in other climatic zones.

3.4. Relationship between Vegetation Dynamics and Climate Factors

3.4.1. Correlation between Vegetation Dynamics with Reversal Trend and Climate Factors

The PCCs between GSN and GS-Tem, GS-Pre, and GS-SR were abbreviated as R_Tem, R_Pre and R_SR, respectively. Before the TP, R_Tem, R_Pre and R_SR were mainly positive correlation, accounting for 58.92%, 68.83% and 56.15%. The regions with R_Tem < 0 were mainly distributed in the northwestern WTSH and the southern MTA. Regions with R_Pre < 0 were mainly distributed in the southern part of plateau climate zone and central WTSH. The regions with R_SR < 0 were mainly concentrated in the central MTA, the western PCSA, the western WTSH and the high altitude region (Figure 7a,c,e). After the TP, the area with R_Tem > 0 and R_SR > 0 decreased from 58.92% to 50.47% and 56.15% to 52.48%, respectively, while the area with R_Pre > 0 increased slightly from 68.83% to 70.25%. Before and after the TP, the areas of R_Tem, R_Pre and R_SR passing the significance test of p < 0.05 increased from 13.26%, 9.73% and 7.38% to 26.48%, 36.52% and 29.13%, respectively (Figure 7b,d,f).

Before and after the TP, the correlation between GSN and the three climate factors had changed significantly (Figure 7a–f). According to the square value of R_Tem, R_Pre and R_SR, the dominant climate factors of vegetation growth were obtained (Figure 7g,h). Before the TP, the temperature-dominated area accounted for 39.34%, followed by precipitation (34.78%). After the TP, vegetation growth was mainly controlled by precipitation (43.70%), and the area controlled by solar radiation (28.83%) was slightly higher than that controlled by temperature (27.47%). The performance of different type vegetation was slightly different. Before the TP, cropland and grassland was mainly affected by temperature (38.57% and 39.11%), followed by precipitation (35.38% and 37.09%), while woodland was mainly affected by temperature and solar radiation (41.49% and 34.88%) (

Table 5. Area percentage of vegetation controlled by different climate factors in different land cover types (%).

Land Cover	Before the Turning Point			After the Turning Point			Monotonic Trend			Non-Significant Trend
	Tem	Pre	SR	Tem	Pre	SR	Tem	Pre	SR	Tem	Pre	SR
Cropland	38.57	35.38	26.05	25.00	48.41	26.59	45.76	39.02	15.22	20.68	57.97	21.35
Woodland	41.49	23.63	34.88	25.27	38.39	36.34	36.12	22.84	41.04	23.77	28.04	48.19
Grassland	39.11	37.09	23.8	28.74	43.40	27.86	52.71	33.90	13.39	25.95	55.15	18.90
ALL	39.34	34.78	25.88	27.47	43.70	28.83	48.04	33.47	18.49	24.99	50.85	24.16

Abbreviations: Tem, Pre, SR represent temperature, precipitation, and solar radiation, respectively.

). After the TP, precipitation became the dominant factor, accounting for 48.41% and 43.40% of cropland and grassland, and the effects of temperature and solar radiation on cropland and grassland were close. The woodland was mainly affected by precipitation and solar radiation (38.39% and 36.34%), and there was little difference in the percentage of their dominant area (Table 5). The climate dominant factor of vegetation growth changed from temperature to precipitation in the northwestern WTSH and the central and southern MTA, while from temperature to solar radiation in PCSH. The correlation between vegetation and climate factors reversed or increased before and after the TP, indicating that the response of vegetation to climate change shifted after the TP.

3.4.2. Correlations between Vegetation Dynamics with Monotonic Trend and Non-Significant Trend and Climate Factors

The GSN with monotonic trend was mainly positively correlated with GS-Tem, GS-Pre and GS-SR with 89.16%, 91.51% and 77.19% the area percentage (Figure 8a,c,e). The area with significant positive correlation between GSN and GS-Tem was slightly higher than that with non-significant positive correlation (Figure 8a). The areas with significant positive correlation between GSN and GS-Pre and GS-SR accounted for 37.01% and 15.77%, respectively (Figure 8c,e). The vegetation growth in the northern and western regions of WTSH, the WTA and plateau climate area was mainly controlled by temperature. Vegetation growth was mainly affected by precipitation in the northwestern part of WTSH and the southern part of MTA. The vegetation growth dominated by solar radiation was mainly distributed in the eastern part of SH, southern of WTSH and CTH (Figure 8g). The significant positive correlation area between NDVI and temperature was higher than that of precipitation and solar radiation. Overall, the vegetation growth with monotonic trend was mainly controlled by temperature, accounting for 48.04% of the study area, followed by precipitation (33.47%) (Figure 8g). The order of main climatic factors affecting cropland and grassland was consistent with the overall vegetation, while woodland was mainly affected by solar radiation (41.04%), followed by temperature (36.12%) (Table 5).

The GSN with non-significant trend was positively correlated with GS-Tem, GS-Pre and GS-SR, accounting for 63.56%, 78.52% and 52.77%, respectively (Figure 8b,d,f). The area with significant positive correlation between GSN and GS-Pre accounted for 31.11%, which was much higher than that between GS-Tem and GS-SR (10.58% and 7.47%, respectively). The regions with significant R_Pre > 0 were mainly distributed in MTSA and MTA. The regions with significant R_SR > 0 were mainly concentrated in the northeastern and southeastern parts of the study area (Figure 8b,d,f). The vegetation with non-significant trend was mainly controlled by precipitation (accounting for 50.85%), which was mainly distributed in arid and semi-arid areas. The vegetation in the plateau climate area was mainly controlled by temperature, while the vegetation dominated by solar radiation was mainly concentrated in the northeastern and southeastern parts of the study area (Figure 8h). Similarly, the impact of precipitation on cropland and grassland was larger than solar radiation and temperature, while solar radiation had the greatest impact on woodland, followed by precipitation (Table 5).

3.5. Identity the Driving Forces in Vegetation Dynamics

Climate change and human activities were conducive to the vegetation growth in most areas, accounting for 76.54% of the total area (Figure 9). A total of 61.50% of the vegetation restoration in the study area was affected by the interaction of climate change and human activities (ICH), of which the average contribution of climate change and human activities to ICH was 48.71% and 51.29%, respectively (Figure 9). ICH-induced vegetation restoration was mainly concentrated in the central, southeastern, southern and western parts of the study area. IHA-induced vegetation restoration was mainly distributed in MTSA, accounting for 5.80% of the study area. However, 9.24% of the area was only affected by ICC, which was distributed in CTH, PCSH and PCSA. The vegetation on 23.46% of the study area had degraded, which was mainly distributed in the northeastern, southwestern and northwestern parts of the study area. DCH was the main driving force of vegetation degradation, accounting for 14.63% of the total area, followed by DHA (5.13%) and DCC (3.70%). The average contribution of climate change and human activities to DCH was 50.79% and 49.21%, respectively (Figure 9).

On the whole, the average contributions of climate change and human activities were 50.33% and 49.67%, with 51.22% and 48.78% for the vegetation restoration area, and 47.53% and 52.57% for the vegetation degradation area, respectively (

Table 6. The contribution rate of climate change (CC) and human activities (HA) to different types of vegetation dynamics in the study area from 1982 to 2015 (%).

Types	Cropland		Woodland		Grassland		ALL Land Cover
	CC	HA	CC	HA	CC	HA	CC	HA
Vegetation restoration	37.10	62.90	57.62	42.38	54.31	45.69	51.22	48.78
Vegetation degradation	56.80	43.20	33.23	66.77	49.72	50.28	47.43	52.57
Overall study area	39.73	60.27	51.04	48.96	53.14	46.86	50.33	49.67

). In general, climate change was the dominant factor in the vegetation restoration area, and human activities were the main driving force in the vegetation degradation area, but the dominant factors of different types of vegetation were significantly different. In the vegetation restoration area, the growth of woodland and grassland was mainly controlled by climatic factors, while the cropland was mainly affected by human activities. In the vegetation degradation area, human activities were the main factors of leading woodland and grassland to degrade, while cropland was mainly controlled by climatic factors. Overall, the contribution of climate change to the dynamics of woodland and grassland was higher than that of human activities, and the contribution of human activities to the dynamics of cropland was much higher than that of climate change (Table 6). The dominant factors of vegetation restoration and degradation had significant spatial heterogeneity (Figure 9a,b). The contribution of human activities was significantly higher than that of climate change in the southeastern part of the study area, while it is generally lower than that of climate change in the western and northwestern parts of the study area.

4. Discussion

This study showed that the trend of NDVI has an increasing rate over time for the EEMD method (Figure 4), which indicated that the process of vegetation was not linear, but nonlinear and variable [10]. Due to its simplicity and convenience, linear trend analysis is still used by many researchers [14,44]. However, linear trend analysis can only obtain the average change rate of vegetation during the study period, and judge the type of change trend [10]. It is easy to underestimate or overestimate the potential risk of vegetation restoration or degradation [12]. The long-term trend obtained by EEMD reveals the change process and trend of vegetation growth over time [23], and has been applied in many studies [6,23,42].

This study showed that the TP of vegetation from browning to greening happened in the early 1990s, which was consistent with the results of Xu et al. [23], while the time of vegetation growth from greening to browning in high altitude areas, such as the southern of the plateau climate zone and the northwestern of the study area, the eastern WTSH, and the CTH was usually later than that in other greening to browning areas (Figure S1). The TP of vegetation has obvious geographical heterogeneity. However, some researchers obtained one or two turning points based on a piecewise linear regression method [45,46] for large-scale regional vegetation change analysis, which will not accurately reflect the real situation of vegetation change. In addition, the vegetation growth trend in Central Asia, Eastern Europe and East Africa [10,12,47,48] had also been reversed on a large scale, which may provide early warning for ecosystem changes under global environmental change [23]; therefore, researchers should pay more attention.

The vegetation in the high altitude area in the northwestern and southwestern parts of the study area experienced a trend transfer process from greening to browning (Figure 5 and Figure 6), which was the same as Liu’s research results [30], but not consistent with the increasing trend considered by other researchers using linear trend analysis [49,50]. From 1982 to 2015, the temperature and solar radiation showed an increasing trend, while the precipitation changed from decreasing to increasing (Figure S2). When the temperature rises, the snow on the mountain will melt, which increases the surface runoff and soil moisture [51,52,53]. In addition, the increase of temperature can improve the root activity of vegetation, improve the photosynthesis rate of vegetation, and contribute to the growth of vegetation [54]. However, the rapid rise in temperature may promote vegetation transpiration, thus accelerating soil water loss [55]. When forced by soil water, vegetation photosynthesis will slow down, which then depresses the growth of vegetation [56]. In addition, the increase in precipitation in high altitude areas means the shortage of sunshine time and the reduction in solar radiation, and the photosynthesis of vegetation will also be suppressed, which is not conducive to the growth of vegetation [57,58]. On the other hand, human activities, such as excessive grassland grazing, are also an important contributor to the vegetation growth trend from greening to browning [59,60]. The vegetation in the central part of the study area, especially the Loess Plateau, had experienced a process from browning trend to greening trend (Figure 5 and Figure 6), which was closely related to ecological projects implemented in China, such as Natural Forest Conservation Program, Grain for Green Program and grazing management, which can promote the significant increase in vegetation coverage [61]. In the early stage, the increase in temperature and solar radiation and the decrease in precipitation triggered the degradation of vegetation in central Inner Mongolia [36]. However, in the later period, the increase in precipitation and the implementation of the Beijing-Tianjin Sand Source Control Program effectively resisted the adverse impact of climate, and promoted the reversal of vegetation growth trend from browning to greening [36,45].

According to different geographical locations and vegetation types, climate change and human activities play different roles in vegetation dynamics [60]. The change rate of cropland is significantly higher than that of woodland and grassland (Figure 4), with cropland generally being distributed in low altitude areas (Figure 1). Therefore, it is highly affected by human activities, such as artificial irrigation, application of chemical fertilizer, and change of planting structure or varieties, which are conducive to vegetation greening [58]. However, the cropland located in the northeastern part of WTSH has a high potential degradation risk (Figure 5 and Figure 6), which may be related to the population increase and the rapid development of urbanization in this area [58]. In the CTH, MTSH and MTSA climate zones, woodland was accompanied by a high potential degradation risk (Table S1). The woodland in CTH was mainly affected by temperature and solar radiation, while in MTSH and MTSA, it was mainly affected by temperature and precipitation [36]. From 1990 to 2006, temperature and solar radiation showed increasing trends, while the precipitation showed a decreasing trend, while from 2006 to 2015, the trend was opposite [36]. The growth of woodland in arid and semi-arid areas requires significantly higher water conditions than that in humid climate areas, so the potential degradation risk of woodland may be caused by adverse factors of climate conditions. In addition, the rapid development of social economy and more urban construction activities may seriously damage the ecological environment, reduce the vegetation coverage, and lead to vegetation degradation [45,62].

The geographical difference of grassland change trend was obvious (Figure 5). Adverse factors of climate change and unreasonable human activities, such as long-term overgrazing, have reduced the self-recovery ability of grassland vegetation and caused the degradation of most natural grasslands to varying degrees [63,64], which seriously affects the sustainable development of grassland animal husbandry and the stability and development of economic and social in pastoral areas. However, favorable climatic conditions and reasonable ecological restoration measures, such as reducing grazing intensity and extending the restoration time of pasture, can effectively improve vegetation productivity and promote vegetation restoration [65]. In addition, human activities can also play a positive role in vegetation restoration under adverse climatic conditions [45]. Therefore, active human intervention and measures can promote the restoration of degraded ecosystems and promote regional sustainable development [66].

The methods used in this study have been widely recognized [12,41,44], but there are still some shortcomings. GIMMS NDVI3g data with 8km spatial resolution cannot reflect more detailed spatial information, and it is difficult to reflect the changes of local vegetation within the pixel scale. Therefore, the NDVI products with higher spatial resolution should be used in future study. Because of the difference in geographical location, the response of vegetation to climate factors has different lag effects [43]. Therefore, in future studies, the lag effects between vegetation and climate factors will be considered in order to improve the accuracy of assessment [41]. In the study, climate factors only selected temperature, precipitation and solar radiation, ignoring the impact of other factors on vegetation change, such as CO₂ [67], nitrogen deposition [68], and soil moisture [69]. The interactive effect of temperature, precipitation and solar radiation is also not considered in this study. In addition, extreme climate events have a great impact on different ecosystems, and the impact on vegetation changes shows regional diversity [70,71]. Therefore, more climate factors will be considered to analyze the impact of climate change on vegetation change, which will help comprehensively understand the relationship between vegetation and climate, and more accurately assess the contribution of human activities to vegetation change. At present, human activities are analyzed as a whole, and how to quantitatively evaluate the impact of different types of human activities on vegetation change is still an important issue [72].

5. Conclusions

This study analyzed the linear and nonlinear trends of vegetation in northern China from 1982 to 2015, assessed the response of different types of vegetation to climate change, and quantified the relative impact of climate change and human activities on vegetation change. The results show that the vegetation in most areas has undergone significant changes during the whole study period. The vegetation growth in 36.64% of the areas had a reverse trend. There were potential risks of vegetation degradation in the southwestern, northwestern and northeastern parts of the study area. The potential degradation risk of woodland was higher than that of cropland and grassland vegetation.

The correlation between vegetation and climate factors was reversed or enhanced before and after the turning point. Temperature and precipitation were the main climatic factors affecting the vegetation with a monotonic trend, and precipitation dominated the vegetation with a non-significant trend. The areas affected by temperature were mainly located in high altitude areas, arid and sub-arid areas were mainly affected by precipitation, and humid and sub-humid areas were mainly affected by solar radiation.

Climate change and human activities had promoted vegetation growth in 76.54% of the study area. Overall, the average contribution rates of climate change and human activities to vegetation change were 50.33% and 49.67%, respectively. The vegetation restoration of grassland and woodland was mainly dominated by climate change, and human activities mainly dominated the vegetation degradation, while the cropland vegetation had the opposite performance. The contribution rate of human activities to the vegetation change in the southeastern and eastern parts of the study area was generally higher than that of climate change, but it was opposite in the higher altitude area, with obvious spatial heterogeneity. The results are helpful to understand the dynamic changes of vegetation in northern China, and provide important reference value for regional ecological protection and restoration.