Detecting Vegetation Variations and Main Drivers over the Agropastoral Ecotone of Northern China through the Ensemble Empirical Mode Decomposition Method

: Vegetation is the major component of the terrestrial ecosystem. Understanding both climate change and anthropogenically induced vegetation variation is essential for ecosystem management. In this study, we used an ensemble empirical mode decomposition (EEMD) method and a linear regression model to investigate spatiotemporal variations in the normalized di ﬀ erence vegetation index (NDVI) over the agropastoral ecotone of northern China (APENC) during the 1982–2015 period. A quantitative approach was proposed based on the residual trend (RESTREND) method to distinguish the e ﬀ ects of climatic (i.e., temperature (TEM), precipitation (PRE), total downward solar radiation (RAD), and near surface wind speed (SWS)) and anthropogenic e ﬀ ects on vegetation variations. The results showed that the NDVI exhibited a signiﬁcant greening trend of 0.002 year − 1 over the entire study period of 1982–2015 and that areas with monotonous greening dominated the entire APENC, occupying 40.97% of the region. A browning trend was also found in the central and northern parts of the APENC. PRE presented the highest spatial correlation with the NDVI and climate factors, suggesting that PRE was the most important factor a ﬀ ecting NDVI changes in the study area. In addition, the RESTREND results indicated that anthropogenic contributions dominated the vegetation variations in the APENC. Therefore, reusing farmland for grass and tree planting made a positive contribution to vegetation restoration, while deforestation, overgrazing, and the reclamation of grasslands were the opposite. In addition, with the continuous implementation of national ecological engineering programs such as the Grain to Green Program, positive human activity contributions to vegetation greening signiﬁcantly increased. These results will support decision- and policy-making in the assessment and rehabilitation of ecosystems in the study region.


Introduction
Agropastoral ecotones (hereafter APEs) refer to transition zones interlaced by croplands and grasslands [1], and they are widely distributed in the arid and semiarid regions of the world [2]. In China, the APEs area spans throughout the whole mainland from northeast to southwest, covering 13.35% of China's land area (1.29 million km 2 ), with residential populations and grain yield accounting for 4.68% (59.30 million) and 7.98% of the total population and grain yield of China, respectively [3]. Therefore, the arid and semiarid agropastoral ecotone of northern China (hereafter APENC) is a typical ecologically transitional zone (largest in area and the longest in spatial scale), which effectively prevents the desert

Study Area Description
The study was conducted across the agropastoral ecotone in northern China (APENC), which lies between 100 • 51 -124 • 35 E and 35 • 3 -47 • 39 N with a general elevation of 72-4,781 m above sea level and has an area of approximately 0.62 million km 2 ( Figure 1). The region is mainly located in the southeastern Inner Mongolia Plateau and the northern Loess Plateau. The APENC is a climate and ecological transition belt, with a semiarid temperate steppe dominated by animal husbandry to the north and a subhumid agricultural plain region to the south that connects the farmland and pastoral areas. The main land use includes grassland, cropland, desert, and forest. The climate of the area is characterized by a transition zone from semiarid to subhumid, with an annual average temperature of approximately 6.0 to 8.5 • C and an annual precipitation ranging from 150 mm in the north to 600 mm in the south. Furthermore, the precipitation is seasonally uneven and mainly occurs during summer in the form of high-intensity rainstorms from July to September, accounting for approximately 60%-70% of the annual precipitation.
Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 22 grasslands, and the Beijing-Tianjin Sandstorm Source Control Project [34,35]. As mentioned above, the state of vegetation in this region is not only related to local ecosystem services and human welfare, but also has an important impact on the environment in northern China. Previous studies have used line regression methods to investigate vegetation variations and related effects of climate and anthropogenic activity [8,36]. However, limitations in the vegetation change trend detection method have hindered our ability to investigate the drivers of vegetation variations [37,38]. Therefore, evaluating vegetation variations using nonlinear methods may help us obtain a deeper understanding of vegetation variations and their correlations with climate change and anthropogenic activities. Thus, the main objectives of this study were to (a) analyze the spatiotemporal vegetation variations in the APENC, (b) identify different vegetation variation trends, (c) separate human-driven from climatically induced vegetation changes, and (d) quantify the positive and negative effects of anthropogenic activity on vegetation variation. In the context of recent pressure from human activities on vegetation, this approach is particularly beneficial for decision-making on land use management and ecological restoration policies.

Study Area Description
The study was conducted across the agropastoral ecotone in northern China (APENC), which lies between 100°51′-124°35′E and 35°3′-47°39′N with a general elevation of 72-4,781 m above sea level and has an area of approximately 0.62 million km 2 ( Figure 1). The region is mainly located in the southeastern Inner Mongolia Plateau and the northern Loess Plateau. The APENC is a climate and ecological transition belt, with a semiarid temperate steppe dominated by animal husbandry to the north and a subhumid agricultural plain region to the south that connects the farmland and pastoral areas. The main land use includes grassland, cropland, desert, and forest. The climate of the area is characterized by a transition zone from semiarid to subhumid, with an annual average temperature of approximately 6.0 to 8.5 °C and an annual precipitation ranging from 150 mm in the north to 600 mm in the south. Furthermore, the precipitation is seasonally uneven and mainly occurs during summer in the form of high-intensity rainstorms from July to September, accounting for approximately 60%-70% of the annual precipitation.

Global Inventory Modeling and Mapping Studies (GIMMS) Data
The third generation of NDVI datasets, which were derived from the Global Inventory Modeling and Mapping Studies (GIMMS NDVI3g) (https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/) and span from 1982 to 2015, were applied to study vegetation variations over the entire APENC. The dataset is a new long-term time series version, with a spatial resolution of 1/12 • and a temporal resolution of 15 days. Yearly NDVI values were calculated through the maximum value composite (MVC) method, which minimizes the influence of clouds, sun angle, water vapor, aerosols, and directional surface reflectance on the NDVI image [39]. To better match the gridded meteorological data, the GIMMS NDVI data were upscaled to 0.1 • × 0.1 • resolution using the bilinear interpolation method in ArcGIS 10.3, and the data sample from 2015 is shown in Figure 2.
and Mapping Studies (GIMMS NDVI3g) (https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/) and span from 1982 to 2015, were applied to study vegetation variations over the entire APENC. The dataset is a new long-term time series version, with a spatial resolution of 1/12° and a temporal resolution of 15 days. Yearly NDVI values were calculated through the maximum value composite (MVC) method, which minimizes the influence of clouds, sun angle, water vapor, aerosols, and directional surface reflectance on the NDVI image [39]. To better match the gridded meteorological data, the GIMMS NDVI data were upscaled to 0.1° × 0.1° resolution using the bilinear interpolation method in ArcGIS 10.3, and the data sample from 2015 is shown in Figure 2.

Climatic Datasets
Climatic datasets, including air temperature, cumulative precipitation, downward shortwave and longwave radiation, and near-surface 10-m wind speed were collected from the China Meteorological Forcing Dataset (CMFD). The entire climate dataset spans from 1979 to 2015, with a temporal resolution of 3 h and a spatial resolution of 0.1° [40]. Previous studies have proven that climatic data from CMFD agree with observation data [41,42]. In this study, the annual precipitation (PRE), average temperature (TEM), near-surface wind speed (SWS), and total downward solar radiation, including shortwave and longwave (RAD), from 1982 to 2015 were used to analyze the correlation between the NDVI and climate factors. The four climate data samples from 2015 are shown in Figure 2.

Climatic Datasets
Climatic datasets, including air temperature, cumulative precipitation, downward shortwave and longwave radiation, and near-surface 10-m wind speed were collected from the China Meteorological Forcing Dataset (CMFD). The entire climate dataset spans from 1979 to 2015, with a temporal resolution of 3 h and a spatial resolution of 0.1 • [40]. Previous studies have proven that climatic data from CMFD agree with observation data [41,42]. In this study, the annual precipitation (PRE), average temperature (TEM), near-surface wind speed (SWS), and total downward solar radiation, including shortwave and longwave (RAD), from 1982 to 2015 were used to analyze the correlation between the NDVI and climate factors. The four climate data samples from 2015 are shown in Figure 2.

Land Use and Land Cover (LUCC) Dataset
The land use and land cover (LUCC) data (1980, 2000, and 2015) used in this study from the Resources and Environmental Data Cloud Platform (http://www.resdc.cn/data.aspx?DATAID=99) [43]. The LUCC dataset was gridded with a 1-km resolution and contains 6 primary land use types, including grasslands, croplands, forests, water bodies, construction lands, and deserts, and 25 secondary land use types are also included. Grasslands are the most widespread land use type in the APENC. To further analyze the land use changes, grassland types with different coverages at the secondary level were used in this study, including high-coverage grasslands, medium-coverage grasslands, and low-coverage grasslands. Finally, eight land cover types, i.e., croplands, forests, water bodies, construction lands, deserts, high-coverage grasslands, medium-coverage grasslands, and low-coverage grasslands, were selected. Notably, considering the time coherence between the different types of datasets in our study, the LUCC data in 1980 were used to substitute the LUCC in 1982 data since such data were not available. Land use and land cover type conversion between 1982 and 2000 and 2000 and 2015 are shown in Figure 3. The land use and land cover (LUCC) data (1980, 2000, and 2015) used in this study from the Resources and Environmental Data Cloud Platform (http://www.resdc.cn/data.aspx?DATAID=99) [43]. The LUCC dataset was gridded with a 1-km resolution and contains 6 primary land use types, including grasslands, croplands, forests, water bodies, construction lands, and deserts, and 25 secondary land use types are also included. Grasslands are the most widespread land use type in the APENC. To further analyze the land use changes, grassland types with different coverages at the secondary level were used in this study, including high-coverage grasslands, medium-coverage grasslands, and low-coverage grasslands. Finally, eight land cover types, i.e., croplands, forests, water bodies, construction lands, deserts, high-coverage grasslands, medium-coverage grasslands, and low-coverage grasslands, were selected. Notably, considering the time coherence between the different types of datasets in our study, the LUCC data in 1980 were used to substitute the LUCC in 1982 data since such data were not available. Land use and land cover type conversion between 1982 and 2000 and 2000 and 2015 are shown in Figure 3. CCR is croplands, FRO is forests, URB is urban construction lands, DES is deserts, HCG is high-coverage grasslands, MCG is medium-coverage grasslands, and LCG is low-coverage grasslands.

Ensemble Empirical Mode Decomposition (EEMD) Method
The EEMD originated from the EMD [44], which is capable of extracting the secular variation trend via the removal of different frequency components from a signal sequence. Overall, the EMD is a dyadic filter capable of separating the time series signal x(t) into different kinds of time series components known as intrinsic mode functions ((IMFs); IMF i = 1 . . . n), and the procedure continues until the residue term R(t) is less than a predetermined value of consequence or the residue term R(t) becomes a monotonic function (Equation (7)) [21]. Then, the signal x(t) is finally expressed as the sum of the (IMF i = 1 . . . n) components plus the final residue term R(t).
However, mode mixing is the obvious defect of the EMD, which implies a single IMF subsisting off other signals from dissimilar scales or a signal of the same scale appearing in divergent IMF components [44]. To overcome the problem of mode mixing in the EMD method, Wu and Huang [25] designed a new noise-assisted signal analysis technique called EEMD, which efficiently improved the mode mixing problem in EMD and preserved the ability of binary filters. The main point of the improved EEMD method is to add white noise to the time series signal x(t) sequence. Through EEMD decomposition, the added white noise can be cancelled out, and the signal x(t) can be decomposed into the IMFs and residue term R(t). The residual term (Rn) is considered to have detached fluctuant components and reserved secular trends, and it represents the intrinsic signal trend. Generally, the Rn is monotonic or contains an extreme value, which is identified as the turning point at the extremum. A flow chart of the proposed method is illustrated in Figure 4, and the decomposition process can be described as follows: Step 1. Add a number of white noise q(t) to original signal o(t) [25]: Step 2. Identify the local extrema of new signal x 1 (t), use cubic spline interpolation to create the upper and lower envelopes, and calculate the mean value e 1 (t) of the upper and lower envelopes [21]: Step 3. Use standard deviation (sd) to determine whether the above process is to be continued or stop sifting [21]: where i is the number of iterations; Step 4. When the sd is smaller than a pregiven value, the above process should be stopped, and the first IMF1 is calculated [21]: Step 5. The rest data r(t) are the differences between IMF1 and signal x 1 (t) [21]: Step 6. Repeat step 1 to step 5 until r n (t) becomes a monotonic function [21]: Remote Sens. 2019, 11, 1860 7 of 23 Step 7. Thus, the number of IMFs and a residual separated from x(t) [21] are obtained: Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 22 In this study, when the EEMD method was applied to analyze the NDVI series, a total of 100 Gaussian white noise series were added to the signal sequence, with an amplitude of 0.2 standard deviations from the signal sequence. Subsequently, the NDVI trend at each pixel was obtained based on the EEMD method, and the accumulated trend was represented as the increment of the NDVI trend compared to the initial time in 1982 (Equation (8)) [45], while the instantaneous trend was regarded as the interannual variation in the NDVI trend (Equation (9)) [45] and the significance test was verified according to Pan's method [9]. In addition, we also calculated the tendencies of TEM, PRE, RAD, and SWS over the entire study period: where t is the end time of the study period (e.g., t = 2015), and t0 is the beginning time of the study period (e.g., t0 = 1982). According to the Rn trend of each NDVI series, the EEMD trends in the NDVI could be classified into five classes. When the NDVI showed a significant monotonic upward or downward trend (P < 0.05), the trend was considered to be monotonic greening or browning. When the trend contained one extreme point and at least one year passed the significance test (P < 0.05), the NDVI trend was classified into greening to browning or browning to greening classes. However, if there was no significant change in the NDVI variation trend, it was classified as no significant change. Four pixels were randomly selected from the NDVI dataset of the APENC as samples of the NDVI variation trend, and they are shown in Figure 5. In this study, when the EEMD method was applied to analyze the NDVI series, a total of 100 Gaussian white noise series were added to the signal sequence, with an amplitude of 0.2 standard deviations from the signal sequence. Subsequently, the NDVI trend at each pixel was obtained based on the EEMD method, and the accumulated trend was represented as the increment of the NDVI trend compared to the initial time in 1982 (Equation (8)) [45], while the instantaneous trend was regarded as the interannual variation in the NDVI trend (Equation (9)) [45] and the significance test was verified according to Pan's method [9]. In addition, we also calculated the tendencies of TEM, PRE, RAD, and SWS over the entire study period: where t is the end time of the study period (e.g., t = 2015), and t 0 is the beginning time of the study period (e.g., t 0 = 1982). According to the Rn trend of each NDVI series, the EEMD trends in the NDVI could be classified into five classes. When the NDVI showed a significant monotonic upward or downward trend (P < 0.05), the trend was considered to be monotonic greening or browning. When the trend contained one extreme point and at least one year passed the significance test (P < 0.05), the NDVI trend was classified into greening to browning or browning to greening classes. However, if there was no significant change in the NDVI variation trend, it was classified as no significant change. Four pixels were randomly selected from the NDVI dataset of the APENC as samples of the NDVI variation trend, and they are shown in Figure 5.

Linear Regression Analysis
The least squares method based on a linear regression was used to investigate the change rate of the NDVI at a pixel scale [46]: where is the NDVI change rate; t is the year order from 1 to n; n is the number of years; and is the NDVI value of year t. A positive (negative) slope indicates a significant increasing (decreasing) vegetation coverage at P < 0.05.

Relationship Analysis between Climate Factors and the NDVI
To evaluate the effect of climatic factors (TEM, PRE, RAD, SWS) on vegetation, we used a partial correlation coefficient to identify the significant relationships between NDVI and climate factors at P < 0.05. In addition, the residual trend (RESTREND) [47] was also applied to separate the climatically induced vegetation changes from the anthropogenically driven vegetation variations. Since the RESTREND method could separate the climate effects from the NDVI trend, the NDVI could be divided into climatically induced (NDVIc) and anthropogenically induced (NDVIa) groups in this study. The specific steps of RESTREND method were as follows [48]: First, a multiple linear regression model was applied to the NDVI and climate variables (TEM, PRE, SAD, SWS) from 1982 to 2015 at each pixel. According to the results obtained from the multiple linear regression model, the NDVIc was calculated. Second, the residual differences between the observed NDVI and predicted NDVIc were calculated as the NDVIa. In addition, the adj-R 2 is generally interpreted as the NDVI variation proportion that can be explained by climate variables. Then, the remaining fraction (i.e., 1adj-R 2 ) was regarded as the anthropogenic effect on NDVI variations. Statistical significance was

Linear Regression Analysis
The least squares method based on a linear regression was used to investigate the change rate of the NDVI at a pixel scale [46]: where Slope NDVI is the NDVI change rate; t is the year order from 1 to n; n is the number of years; and NDVI t is the NDVI value of year t. A positive (negative) slope indicates a significant increasing (decreasing) vegetation coverage at P < 0.05.

Relationship Analysis between Climate Factors and the NDVI
To evaluate the effect of climatic factors (TEM, PRE, RAD, SWS) on vegetation, we used a partial correlation coefficient to identify the significant relationships between NDVI and climate factors at P < 0.05. In addition, the residual trend (RESTREND) [47] was also applied to separate the climatically induced vegetation changes from the anthropogenically driven vegetation variations. Since the RESTREND method could separate the climate effects from the NDVI trend, the NDVI could be divided into climatically induced (NDVIc) and anthropogenically induced (NDVIa) groups in this study. The specific steps of RESTREND method were as follows [48]: First, a multiple linear regression model was applied to the NDVI and climate variables (TEM, PRE, SAD, SWS) from 1982 to 2015 at each pixel. According to the results obtained from the multiple linear regression model, the NDVIc was calculated. Second, the residual differences between the observed NDVI and predicted NDVIc were calculated as the NDVIa. In addition, the adj-R 2 is generally interpreted as the NDVI variation proportion that can be explained by climate variables. Then, the remaining fraction (i.e., 1 -adj-R 2 ) was regarded as the anthropogenic effect on NDVI variations. Statistical significance was tested using a t-test, and only grids with significant (P <0.05) temporal trends were used for further analysis.

Relative Contribution Rates of Anthropogenic Factors to the NDVI
To distinguish and quantify the anthropogenic effects on the NDVI, we proposed a quantitative approach as an extension of Liu's method [49]. Based on the above testing and calculations, the grids with 95% significance values were used in the following processing. The contribution of anthropogenic factors was calculated by Equation (11): where C_NDVI a refers to the anthropogenic contribution, NDVI t n a refer to every year except the first year NDVI of the anthropogenic effect and NDVI t 0 a refer to the first year NDVI of the anthropogenic effect, respectively. When C_NDVI a > 0, this indicates that anthropogenic activities have positive effects on vegetation (e.g., conservation and restoration measures), and we use C_NDVI ap to express the positive anthropocentricity. In contrast, when C_NDVI a < 0, anthropocentric factors have negative effects on vegetation (e.g., overgrazing, deforestation), and we use C_NDVI an to express the negative anthropocentricity. Therefore, the total contributions of positive NDVIa (SUM_C_NDVI ap ) and negative NDVIa (SUM_C_NDVI an ) were calculated by Equation (12) and Equation (13): In this research, since large-scale revegetation activities were implemented in this area around 2000, the whole research period was divided into two sections to explore the impacts of anthropogenic activities on vegetation variation from 1982 to 2000 and from 2000 to 2015 (on a pixel scale). Therefore, during 1982-2000, the t 0 and t n referred to 1982 and 2000, while during 2000-2015, the t 0 and t n referred to 2000 and 2015, respectively. Then, the contribution rates of anthropogenic activities and land use were further superimposed to calculate the contributions of different land use changes to the NDVI. The whole workflow is shown in Figure 6.
Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 22 tested using a t-test, and only grids with significant (P <0.05) temporal trends were used for further analysis.

Relative Contribution Rates of Anthropogenic Factors to the NDVI
To distinguish and quantify the anthropogenic effects on the NDVI, we proposed a quantitative approach as an extension of Liu's method [49]. Based on the above testing and calculations, the grids with 95% significance values were used in the following processing. The contribution of anthropogenic factors was calculated by Equation (11): where _ refers to the anthropogenic contribution, refer to every year except the first year NDVI of the anthropogenic effect and refer to the first year NDVI of the anthropogenic effect, respectively. When _ > 0, this indicates that anthropogenic activities have positive effects on vegetation (e.g., conservation and restoration measures), and we use _ to express the positive anthropocentricity. In contrast, when _ < 0, anthropocentric factors have negative effects on vegetation (e.g., overgrazing, deforestation), and we use _ to express the negative anthropocentricity. Therefore, the total contributions of positive NDVIa ( _ _ ) and negative NDVIa ( _ _ ) were calculated by Equation (12) and Equation (13): In this research, since large-scale revegetation activities were implemented in this area around 2000, the whole research period was divided into two sections to explore the impacts of anthropogenic activities on vegetation variation from 1982 to 2000 and from 2000 to 2015 (on a pixel scale). Therefore, during 1982-2000, the t0 and tn referred to 1982 and 2000, while during 2000-2015, the t0 and tn referred to 2000 and 2015, respectively. Then, the contribution rates of anthropogenic activities and land use were further superimposed to calculate the contributions of different land use changes to the NDVI. The whole workflow is shown in Figure 6.
Remote Sens. 2019, 11, x FOR PEER REVIEW 10 of 22 Figure 6. Workflow of the NDVI analysis methods and processes.

Temporal and Spatial Trends in the NDVI in the APENC
The linear and EEMD variation trends in the NDVI in the APENC are presented in Figure 7. The linear trend suggested that the annual mean NDVI of the APENC increased at a rate of 0.001 year -1 (P < 0.01) during the 1982-2015 period. Compared to the NDVI before the 2000s (1982-2000), the NDVI trend rate increased by 0.0013 year -1 during the 2000s (Figure 7a). Similarly, the EEMD trend presented a significant increasing trend (slope = 0.002 year -1 , P < 0.05) over the entire study period (Figure 7b). Nonetheless, the instantaneous greening rate first decreased and then increased, with a corresponding rate of decline from 0.002 year -1 to 0.001 year -1 during 1982-2000 than increased to 0.003 year -1 in 2015 (Figure 7c). In general, the NDVI in the APENC continuously increased from 1982 to 2015, and the greening rate increased significantly after 2000.  Figure 8a presents the spatial patterns of the NDVI variation trend in the APENC from 1982 to 2015. The linear trend showed that approximately 49.03% of the areas displayed a significant increase in the NDVI (slope = 0.003 year -1 , P < 0.05), and these areas were primarily located in the southwest and (partially) northeast of the study region, especially in the Shaanxi, Ningxia, and Gansu provinces. However, nearly 9.48% of the region had a significant vegetation browning tendency (slope = -0.002 year -1 , P < 0.05), which was mainly concentrated in the northeastern APENC. The remaining areas (accounting for 41.49% of the total region), distributed throughout the northcentral part of the APENC, showed no remarkable changes. Over the entire time series, approximately 96.24% of the regional NDVI that showed significant trends was detected through EEMD (Figure 8b), including 74.60% of the area experiencing greening (slope = 0.003 year -1 , P < 0.05) and 21.64% of the area experiencing browning (slope = 0.001 year -1 , P < 0.05). The spatial distribution of the EEMD trend was generally unanimous with a linear trend, except for the northcentral part of the research region, where the magnitude of the vegetation change rate was relatively low (approximately ±0.002).
In addition, the spatial trend classification of EEMD is shown in Figure 8c. Considering the spatial pattern of the EEMD trends, greening to greening was the main type, accounting for 40.97%  Figure 8a presents the spatial patterns of the NDVI variation trend in the APENC from 1982 to 2015. The linear trend showed that approximately 49.03% of the areas displayed a significant increase in the NDVI (slope = 0.003 year −1 , P < 0.05), and these areas were primarily located in the southwest and (partially) northeast of the study region, especially in the Shaanxi, Ningxia, and Gansu provinces. However, nearly 9.48% of the region had a significant vegetation browning tendency (slope = −0.002 year −1 , P < 0.05), which was mainly concentrated in the northeastern APENC. The remaining areas (accounting for 41.49% of the total region), distributed throughout the northcentral part of the APENC, showed no remarkable changes. Over the entire time series, approximately 96.24% of the regional NDVI that showed significant trends was detected through EEMD (Figure 8b), including 74.60% of the area experiencing greening (slope = 0.003 year −1 , P < 0.05) and 21.64% of the area experiencing browning (slope = 0.001 year −1 , P < 0.05). The spatial distribution of the EEMD trend was generally unanimous with a linear trend, except for the northcentral part of the research region, where the magnitude of the vegetation change rate was relatively low (approximately ±0.002).
In addition, the spatial trend classification of EEMD is shown in Figure 8c. Considering the spatial pattern of the EEMD trends, greening to greening was the main type, accounting for 40.97% of the vegetated area. In comparison, areas that kept browning were smaller, occupying only approximately 4.49% of the vegetated area, with a mean annual browning rate of 0.0024 year −1 (P < 0.05). In addition to monotonic greening or browning, more than 49.96% of the vegetated area showed a reversed vegetation growth tendency. Roughly 26.78% and 23.99% of the vegetated area were greening to browning and browning to greening, respectively. Spatially, the greening to browning was concentrated in the northcentral part, and the browning to greening was scattered throughout the study region. of the vegetated area. In comparison, areas that kept browning were smaller, occupying only approximately 4.49% of the vegetated area, with a mean annual browning rate of 0.0024 year -1 (P < 0.05). In addition to monotonic greening or browning, more than 49.96% of the vegetated area showed a reversed vegetation growth tendency. Roughly 26.78% and 23.99% of the vegetated area were greening to browning and browning to greening, respectively. Spatially, the greening to browning was concentrated in the northcentral part, and the browning to greening was scattered throughout the study region.  Figure 9 displays a comparison between the EEMD trend distribution and the NDVI variation rate before and after 2000. The results showed that the NDVI variation rate after 2000 (slope = 0.002 year -1 , P < 0.05) was significantly higher than before 2000 (slope = 0.001 year -1 , P < 0.05), and the areas experiencing greening to greening enlarged remarkably. However, a more noteworthy phenomenon was that browning to browning was constantly expanding around the northcentral part of the APENC, suggesting an increasing degradation rate. Overall, there were significant differences between the western and eastern parts of the study area in terms of vegetation change. That is, the grassland variability in the western part was greater than that in the eastern APENC.  Figure 9 displays a comparison between the EEMD trend distribution and the NDVI variation rate before and after 2000. The results showed that the NDVI variation rate after 2000 (slope = 0.002 year −1 , P < 0.05) was significantly higher than before 2000 (slope = 0.001 year −1 , P < 0.05), and the areas experiencing greening to greening enlarged remarkably. However, a more noteworthy phenomenon was that browning to browning was constantly expanding around the northcentral part of the APENC, suggesting an increasing degradation rate. Overall, there were significant differences between the western and eastern parts of the study area in terms of vegetation change. That is, the grassland variability in the western part was greater than that in the eastern APENC. Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 22 Figure 9. Spatial distribution of the EEMD trend classification (a,c) and variation rate (b,d) before and after 2000. Only pixels with P < 0.05 are presented.

Relationships between Vegetation Dynamics and Climate Change
The temporal and spatial linear and EEMD trends in climatic factors are shown in Figure 10 and Figure 11. TEM, PRE, and RAD increased by 0.05 °C/year (P < 0.05), 3.54 mm/year (P < 0.05), and 1.43 MJ/m 2 /year (P < 0.05), respectively, during 1982-2015. Both TEM and RAD increased rapidly during 1982-2000, while after 2000, TEM was maintained at a high level and RAD declined slightly. PRE generated larger interannual fluctuations, with mean annual precipitation changes between 283.67 and 539.60 mm. The EEMD trend showed that PRE fluctuated during the entire study period, with variations that first decreased and then increased between 1982 and 2000, and a continuously increasing pattern was observed after 2000. Compared to the general upward trend in the other three meteorological factors, SWS maintained a monotonic decline during 1982-2015. Spatially, the variation trend in these climatic factors revealed an obvious spatial heterogeneity. TEM, RAD, and PRE varied with a mostly obvious zonal distribution in space. A higher warming trend occurred in the southwest, and apparent decreasing precipitation was discovered in the northcentral part of the study region. Over the whole period, SWS decreased by -0.02 m/year (P < 0.05), and a dramatic decrease was detected in the northcentral and northeastern parts of the APENC.

Relationships between Vegetation Dynamics and Climate Change
The temporal and spatial linear and EEMD trends in climatic factors are shown in Figures 10 and 11. TEM, PRE, and RAD increased by 0.05 • C/year (P < 0.05), 3.54 mm/year (P < 0.05), and 1.43 MJ/m 2 /year (P < 0.05), respectively, during 1982-2015. Both TEM and RAD increased rapidly during 1982-2000, while after 2000, TEM was maintained at a high level and RAD declined slightly. PRE generated larger interannual fluctuations, with mean annual precipitation changes between 283.67 and 539.60 mm. The EEMD trend showed that PRE fluctuated during the entire study period, with variations that first decreased and then increased between 1982 and 2000, and a continuously increasing pattern was observed after 2000. Compared to the general upward trend in the other three meteorological factors, SWS maintained a monotonic decline during 1982-2015. Spatially, the variation trend in these climatic factors revealed an obvious spatial heterogeneity. TEM, RAD, and PRE varied with a mostly obvious zonal distribution in space. A higher warming trend occurred in the southwest, and apparent decreasing precipitation was discovered in the northcentral part of the study region. Over the whole period, SWS decreased by −0.02 m/year (P < 0.05), and a dramatic decrease was detected in the northcentral and northeastern parts of the APENC. Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 22   We investigated the effects of climate factors (TEM, PRE, RAD, and SWS) on NDVI variations over the entire period, and the spatial distributions of the partial correlation between the vegetation and climatic factors are shown in Figure 12. Throughout the study period, over 47.79% of the area showed that PRE had a positive effect (P < 0.05) on vegetation greening, and the mean partial correlation coefficient was more than 0.46. These areas were largely located in nonforest and nondesert areas. The positive correlation for these areas suggested that PRE was a major influential factor limiting vegetation growth. In addition, the SWS was the second most important climatic factor affecting vegetation: approximately 23.43% of the study area showed a negative correlation, mainly distributed in the northeast part of the study region. The impacts of TEM on vegetation displayed distinct north-south differences, with significant positive correlations occurring in the southwestern part (9.14%), which is dominated by cultivated farmland. A negative correlation was mainly concentrated in the northeastern part (4.01%), where grassland is the main land use type. Compared to PRE, SWS, and TEM, the effect of RAD on NDVI was quite small, with only 5.63% and 2.05% of the total area showing positive and negative correlations, respectively.
To further reveal the effect of climate on vegetation variations, we compared the partial correlation before and after 2000. For PRE, the positive effect of precipitation after 2000 was more evident than before 2000, and it increased significantly in the northeastern part of the study area. With a continuous increase in TEM, negative temperature effects on vegetation gradually appeared, especially in the central and northern parts of the APENC. The negative correlation range of temperature constantly expanded, although the correlation was not significant. We compared the effects of the SWS on vegetation before and after 2000: as the SWS decreased, the positive effects of the SWS on vegetation increased. We investigated the effects of climate factors (TEM, PRE, RAD, and SWS) on NDVI variations over the entire period, and the spatial distributions of the partial correlation between the vegetation and climatic factors are shown in Figure 12. Throughout the study period, over 47.79% of the area showed that PRE had a positive effect (P < 0.05) on vegetation greening, and the mean partial correlation coefficient was more than 0.46. These areas were largely located in nonforest and nondesert areas. The positive correlation for these areas suggested that PRE was a major influential factor limiting vegetation growth. In addition, the SWS was the second most important climatic factor affecting vegetation: approximately 23.43% of the study area showed a negative correlation, mainly distributed in the northeast part of the study region. The impacts of TEM on vegetation displayed distinct north-south differences, with significant positive correlations occurring in the southwestern part (9.14%), which is dominated by cultivated farmland. A negative correlation was mainly concentrated in the northeastern part (4.01%), where grassland is the main land use type. Compared to PRE, SWS, and TEM, the effect of RAD on NDVI was quite small, with only 5.63% and 2.05% of the total area showing positive and negative correlations, respectively.
To further reveal the effect of climate on vegetation variations, we compared the partial correlation before and after 2000. For PRE, the positive effect of precipitation after 2000 was more evident than before 2000, and it increased significantly in the northeastern part of the study area. With a continuous increase in TEM, negative temperature effects on vegetation gradually appeared, especially in the central and northern parts of the APENC. The negative correlation range of temperature constantly expanded, although the correlation was not significant. We compared the effects of the SWS on vegetation before and after 2000: as the SWS decreased, the positive effects of the SWS on vegetation increased.

Relationships between Vegetation Dynamics and Anthropogenic Activity
In addition to climate, anthropogenic activity is also highly important in terms of affecting vegetation variation. Based on the results of a multivariate regression analysis, the effects of climatic and anthropogenic factors on vegetation during 1982-2015 were distinguished over the entire study region. The numerical results indicated that approximately 53.85% of the APENC passed the significance test (P < 0.05), where the areas that experienced the strongest influences from climate and anthropogenic activities were distributed in the southwestern and northeastern parts of the APENC, respectively. In general, human activity contributed approximately 63.85% to vegetation variation (Figure 13).

Relationships between Vegetation Dynamics and Anthropogenic Activity
In addition to climate, anthropogenic activity is also highly important in terms of affecting vegetation variation. Based on the results of a multivariate regression analysis, the effects of climatic and anthropogenic factors on vegetation during 1982-2015 were distinguished over the entire study region. The numerical results indicated that approximately 53.85% of the APENC passed the significance test (P < 0.05), where the areas that experienced the strongest influences from climate and anthropogenic activities were distributed in the southwestern and northeastern parts of the APENC, respectively. In general, human activity contributed approximately 63.85% to vegetation variation ( Figure 13). We also compared the human activity contribution rate before and after 2000 ( Figure 14). The positive human activity contribution rate increased significantly after 2000, with the average contribution rate increasing from 44.16% to 54.53%, especially in the Loess Plateau located in the southwest of the APENC. On the contrary, the negative human activity contribution rate decreased from 54.54% to 45.46%. Compared to the southwestern part of the study area, the positive human activity contribution rate in the northeastern APENC was much lower, while the negative contribution of human activities increased in some areas. We also compared the human activity contribution rate before and after 2000 ( Figure 14). The positive human activity contribution rate increased significantly after 2000, with the average contribution rate increasing from 44.16% to 54.53%, especially in the Loess Plateau located in the southwest of the APENC. On the contrary, the negative human activity contribution rate decreased from 54.54% to 45.46%. Compared to the southwestern part of the study area, the positive human activity contribution rate in the northeastern APENC was much lower, while the negative contribution of human activities increased in some areas. Remote Sens. 2019, 11, x FOR PEER REVIEW 16 of 22  Table 1. During 1982-2000, the negative effects of anthropogenic activities occupied a significant proportion. In the areas without land use change, persistent cropland and persistent lowand medium-coverage grassland were negatively affected by human activities, and the negative contributions were near 60% of the total anthropogenic activities. In the region that underwent land use changes, the major land use change was the conversion of grasslands (including high-, medium-, and low-coverage) to croplands, which covered an area of approximately 7.0 × 10 3 km 2 . The land use change from low-coverage grasslands to croplands was seriously influenced by human activities. However, after 2000, the negative effects from anthropogenic activities decreased, and the positive effects increased, accounting for 54.22% and 45.78% of the total anthropogenic disturbance, respectively. In the areas with no land use change, croplands, water bodies, construction lands, and medium-and low-coverage grasslands were positively affected by human activities, and the positive contribution of anthropogenic activities was more than 55% of the total anthropogenic activities. In contrast, in the areas with land use changes, the conversions between grasslands and croplands were still the main land use change type: Approximately 1.1 × 10 3 km 2 of the croplands area was transformed to grasslands, and 1.2 × 10 3 km 2 of the grassland area was converted to croplands, which accounted for nearly 16% of the total area during 1980-2000. Of these changes, the contributions of positive human activities caused by land use conversion from croplands to grasslands and mediumcoverage grasslands to croplands were more than 57% and 54%, respectively.   Table 1. During 1982-2000, the negative effects of anthropogenic activities occupied a significant proportion. In the areas without land use change, persistent cropland and persistent low-and medium-coverage grassland were negatively affected by human activities, and the negative contributions were near 60% of the total anthropogenic activities. In the region that underwent land use changes, the major land use change was the conversion of grasslands (including high-, medium-, and low-coverage) to croplands, which covered an area of approximately 7.0 × 10 3 km 2 . The land use change from low-coverage grasslands to croplands was seriously influenced by human activities. However, after 2000, the negative effects from anthropogenic activities decreased, and the positive effects increased, accounting for 54.22% and 45.78% of the total anthropogenic disturbance, respectively. In the areas with no land use change, croplands, water bodies, construction lands, and medium-and low-coverage grasslands were positively affected by human activities, and the positive contribution of anthropogenic activities was more than 55% of the total anthropogenic activities. In contrast, in the areas with land use changes, the conversions between grasslands and croplands were still the main land use change type: Approximately 1.1 × 10 3 km 2 of the croplands area was transformed to grasslands, and 1.2 × 10 3 km 2 of the grassland area was converted to croplands, which accounted for nearly 16% of the total area during 1980-2000. Of these changes, the contributions of positive human activities caused by land use conversion from croplands to grasslands and medium-coverage grasslands to croplands were more than 57% and 54%, respectively.

Advantages of the EEMD Method in Decomposition of the Vegetation Variation Trend
EEMD is among the advanced time-frequency analysis methods suitable for nonlinear and nonstationary signals in many fields, and EEMD has been introduced into vegetation dynamic research [20]. In this study, we used both linear regression and nonlinear EEMD methods to reveal the NDVI variations in the APENC from 1982 to 2015. Both methods showed that the NDVI generally increased during this period, which was consistent with the results of Liu [8]. Comparing the two methods, linear regression can directly describe the overall trend in vegetation dynamics. However, regarding nonlinear relationships between variables, linear regressions are unable to flexibly capture complex change patterns, such as those in some nonmonotonic variation areas. In contrast, the EEMD method can extract the intrinsic time trend from the long-time signal sequence and remove noise and disturbances on an interannual scale [9]. Thus, the EEMD method reflected the browning trend hidden under the greening trend in NDVI, showing that the interannual greening rate of vegetation slowed and that the area of greening to browning increased continuously from 1982 to 2000, while from 2000 to 2015, the greening trend increased and the browning trend gradually decreased. In addition, the EEMD method could break down an NDVI series into its characteristic times and apply a turning point to the temporal trend classification of the NDVI. Even in the northeastern part of the APENC, where the vegetation variation rate was small (roughly from −0.02 to 0.02 year −1 ) (Figure 8), the EEMD method could still capture the vegetation change trend. The classification of this variation trend is helpful for understanding vegetation dynamics more intuitively.

The NDVI Response to Climate Change
Since the APENC is located on the edge of the Asian monsoon and is subjected to both the Asian Monsoon and the Western Current [50], this area is highly sensitive to climatic oscillations [51]. During the last 34 years, dramatic climate changes have occurred over the entire APENC, showing a remarkable upward tendency in the TEM, PRE, and RAD and a downward tendency in the SWS (Figure 11). The climate in northern China has changed from warm-dry to warm-wet since the late 1980s [52]. These climatic changes have observably altered plant growth [53], and the PRE influence on vegetation variation was greater than that of TEM, RAD, and SWS. Previous studies have reported that precipitation is the main limiting climatic factor affecting vegetation growth in arid and semiarid regions [54,55]. An increased PRE tendency likely results in sufficient soil moisture supplementation, but the trends in TEM and solar radiation can induce much more plant transpiration and soil evaporation [8]. Therefore, in the areas dominated by grasslands in the northeast of the APENC, for example, the TEM increase exhibited a negative effect on vegetation, especially during 1982-2000 when the TEM and RAD increased significantly, but the PRE presented a decreasing trend. As a result, a large amount of vegetation in the northeastern APENC was degraded, with most of the vegetation experiencing greening to browning during 1982-2000, with continued browning after 2000 ( Figure 12). This degradation phenomenon did not weaken until the precipitation gradually increased after 2005, and this reversal phenomenon was also found in the Mongolian Plateau [56]. In contrast, in the southwestern APENC, which is dominated by arable land, the warming trend extended the grain-filling duration and thereby increased crop production [57]. On the other hand, the irrigation of arable land could also have reduced the vapor pressure deficit. In addition to precipitation, SWS was the second crucial climatic factor influencing vegetation variation, since the wind speed was negatively correlated with soil moisture content, soil texture, and nutrient content [58][59][60]. A similar result indicated that the SWS has shown a downward trend in recent decades [61]. Therefore, the decreased SWS in recent decades had a positive effect on vegetation growth from a soil moisture perspective.

NDVI Response to Anthropogenic Activities
Vegetation variation is also significantly impacted by anthropogenic activities. Accurate quantification of the impacts of anthropogenic disturbances on vegetation is beneficial to comprehensively understand vegetation variations. The residual trend method (RESTREND) has been widely used to separate anthropogenic activities from climatic effects on vegetation by developing an NDVI climatic model. Li et al. [34] and Liu et al. [49] constructed models to assess the relative contribution rates of unchanged and changed land uses to vegetation variation. Based on previous research, we improved the method to identify the effects of land use change on vegetation to further distinguish between the positive and negative effects of climatic and anthropogenic disturbances. The results suggested that anthropogenic contributions dominated the vegetation variations, which was consistent with the results of Liu et al. [8]. Before 2000, a large amount of grassland was converted to cropland, and the high-coverage grassland was degraded to low-coverage grassland, which indicated that human activities played a negative role in vegetation change. Since the economic reforms in China in 1978, croplands and livestock owned by the collectives or human communes have been redistributed to individual households [62,63]. According to economic statistics, the number of livestock increased significantly during 1982-2000 [64], which contemporaneously contributed to an obvious decrease in the NDVI magnitude. Therefore, extensive cultivation and overgrazing were important causes of vegetation degradation during 1982-2000.
After 2000, a series of land use changes were implemented in the degraded grasslands, including land use conversions from croplands and desert to grasslands and from low-coverage grassland to high-coverage grassland. The positive human contribution rate significantly increased, especially in the Loess Plateau located in the southwestern APENC. These land use changes were attributed to the implementation of ecological restoration strategies (i.e., the Grain to Green Program and the Returning Grazing Land to Grassland Policy), and land rehabilitation has been obviously increasing since the 2000s [65]. In 2003, "grassland transfer" policies were proposed in agricultural and pastoral areas in many provinces, including Inner Mongolia, Shaanxi, Ningxia, Gansu, and Xinjiang [42]. Based on these policies, grasslands gradually underwent rehabilitation: Both high-coverage and medium-coverage grasslands tended to increase [66], and croplands and desert were returned to grasslands and wetlands gradually [67]. Therefore, the negative human contribution rate decreased in the northeast of the APENC. However, a significant cropland expansion was also detected in the northeastern part of the study area and was mainly converted from grasslands. Wei et al. [68] indicated that the abolition of an agricultural tax and the implementation of a farmland protection system by the Chinese government were crucial reasons for cropland areas rebounding.

Limitations of the Study
In this study, we first used the RESTREND method to separate the vegetation changes induced by climatic or anthropogenic factors and further distinguished the positive and negative contributions of anthropogenic activities. The results revealed some differences in the effects of climatic factors and anthropogenic activities on vegetation change. Although the RESTREND method has been extensively applied to study vegetation change, the method still contains some drawbacks. On the one hand, in the establishment of the regression equation of climatic factors and the NDVI, the proportions of pixels with statistical significance showed an important relationship with the time period of the regression equation [69]. To solve this problem, in our study, we attempted to use a longer time scale to reduce the method's biases. On the other hand, the selection of climatic factors had an important influence on the inverse mode of the climate-driven NDVI. In this paper, we chose four major climatic factors (PRE, TEM, RAD and SWS), and of course, other climatic factors may also have had an impact on the vegetation change in this region. Therefore, in future studies, more climatic factors need to be applied to analyze climate effects on NDVI variations, which may increase the proportion of the NDVI explained by climatic factors via linear regression. Considering additional climatic factors can also help to comprehensively understand the relationship between vegetation and climate and more accurately assess the contribution of human activity to vegetation change.

Conclusions
In this study, we used the EEMD method to assess the spatiotemporal dynamics of vegetation variation in the agropastoral ecotone of northern China during 1982-2015. Subsequently, we also distinguished and quantified the relative effects of climatic and anthropogenic factors on vegetation variation. The results showed that most areas of the APENC experienced a monotonous greening trend throughout the entire study period. However, under the background of the overall greening trend, there was a notable trend of vegetation browning in the northern and central parts of the study area. The NDVI and precipitation presented a remarkable spatiotemporal correlation, suggesting that precipitation was the main climatic driver over the whole area. Additionally, a decrease in precipitation and an increase in temperature were two important reasons for vegetation degradation in the northcentral part of the APENC. Subsequently, the RESTREND results indicated that anthropogenic activity was the major driving factor (63.85%) of vegetation variation, and the positive human contribution rate significantly increased after 2000, especially in the southwest of the APENC.