Identification of Natural and Anthropogenic Drivers of Vegetation Change in the Beijing-Tianjin-Hebei Megacity Region

Abstract

Identifying the natural and anthropogenic mechanisms of vegetation changes is the basis for adapting to climate change and optimizing human activities. The Beijing-Tianjin-Hebei megacity region, which is characterized by significant geomorphic gradients, was chosen as the case study area. The ordinary least squares (OLS) method was used to calculate the NDVI trends and related factors from 2000 to 2015. A geographic weighted regression (GWR) model of NDVI trends was constructed using 14 elements of seven categories. Combined with the GWR calculation results, the mechanisms of the effects of explanatory variables on NDVI changes were analyzed. The findings suggest that the overall vegetation displayed an increasing trend from 2000 to 2015, with an NDVI increase of ca. 0.005/year. Additionally, the NDVI fluctuations in individual years were closely related to precipitation and temperature anomalies. The spatial pattern of the NDVI change was highly consistent with the gradients of geomorphology, climate, and human activities, which have a tendency to gradually change from northwest to southeast. The dominant climate-driven area accounted for only 5.98% of the total study area. The vegetation improvement areas were regionally concentrated and had various driving factors, and vegetation degradation exhibited strong spatial heterogeneity. The vegetation degradation was mainly caused by human activities. Natural vegetation was improved because of natural factors and reductions in human activities. Moreover, cropland vegetation as well as urban and built-up area improvements were related to increased human actions and decreased natural effects. This study can assist in ecological restoration planning and ecological engineering implementation in the study area.

1. Introduction

Under the combined effects of global climate change and regional socioeconomic activities, vegetation at different spatial scales is undergoing complex change processes [1,2]. Thus, rationally distinguishing the natural and anthropogenic drivers of vegetation change has become the basis for ecological restoration. The normalized difference vegetation index (NDVI) is an important parameter of the vegetation environment [3,4], and its spatial differentiation and dynamics have become important indicators for ecological environment monitoring [5]. Therefore, many researchers have applied different types and multi-temporal NDVI products [6,7,8] to study the patterns and driving forces of vegetation changes.

At a global scale, vegetation differentiation and trends are affected by various factors, such as climate [9,10], landforms [11,12], land cover types [13,14], and human activities [14,15,16,17]. For example, de Jong et al. [18] believed that global terrestrial vegetation activities had generally declined with climate change from 1982 to 2008. Strong ENSO events and volcanic activity made significant effects on dramatic changes in vegetation. Zhang et al. [19] found that global vegetation showed a remarkable greening trend from 2000 to 2015 based on MODIS-C6 vegetation index product, and stronger association between the climate and vegetation in the boreal and arid regions than other areas. Chen et al. [14] pointed out that global MODIS data revealed increasing leaf area of vegetation from 2000 to 2017, which is caused by human land-use management and climate change, CO₂ fertilization, nitrogen deposition, and recovery from natural disturbances. Liu et al. [15] found that global NDVI increased from 1982 to 2012. However, it is unclear how natural and human factors together affect NDVI trends and how their effects evolve over time.

At the regional scale, especially in urban agglomerations, natural and anthropogenic factors are strongly coupled [20]. Researchers face the interaction of them [21], especially multidimensional and multivariant human activities [22], because various factors have different effects on NDVI differentiation and trends, such as convergence, divergence, nonlinearity, threshold, and feedback [23,24,25]. For example, Barrera and Henríquez [26] studied three urban agglomerations in Chile located in different climates found that urban expansion led to reduced vegetation from 1986 to 2015, and also promoted the improvement of vegetation inside the cities [26]. Zewdie et al. [27] monitored the NDVI dynamics in Northwestern Ethiopia from 2000 to 2014, which found that climate variables were not the main explaining factors, and the decline in vegetation was affected by the increasing pressure of human activities. Peng et al. [28] found the NDVI generally increased in Eastern China from 1999 to 2008 with complex differences driven by precipitation and temperature, and the impact of anthropogenic activities on vegetation dynamics had accumulative effects and a phase effect. Hu et al. [29] found that the NDVI overall improved from 2000 to 2016, which was influenced by climate change, topographic factors, afforestation, resident lifestyles, and economic development policies in the Pearl River Delta of China.

Previous research has revealed the differentiation pattern, change dynamics, and driving mechanism from various angles. However, there is no consensus on the types, strengths, spatial distributions, and time series evolution of natural and anthropogenic drivers of vegetation changes at the global or regional scale, and impacts such as geomorphic gradients and land cover types have not been fully considered. Can we solve these problems to a certain extent through reasonable geographic division and matching data selection and regional scale method application? We chose the Beijing-Tianjin-Hebei megacity region with significant geomorphic gradients as the research area, which is experiencing high-speed economic development and facing complex ecological problems. This study was performed to (i) detect the overall trend of the NDVI changes from 2000 to 2015, (ii) reveal the spatial differentiation characteristics of the NDVI trends, and (iii) distinguish the natural and anthropogenic drivers of NDVI trends in different land cover types and different geomorphic units.

2. Materials and Methods

2.1. Study Area

The Beijing-Tianjin-Hebei megacity region (abbreviated as BTH) is China’s political and cultural center and an important core area for the economy of Northern China (Figure 1a). The geographical scope (GCS_WGS_1984) is 113.457–119.8526°E, 36.0431–42.6216°N, which covers an area of approximately 218,000 km². The region belongs to the continental semi-humid and semiarid monsoon climate zone, where the annual mean precipitation was 534 mm, and the average temperature was 9.87 ℃ from 1980–2015, according to the data from the Resources and Environmental Data Cloud Platform (http://www.resdc.cn). The types of vegetation in BTH include natural vegetation and crops, of which natural vegetation includes warm temperate deciduous broad-leaved forests and temperate grasslands, and crops mainly include wheat, corn, and cotton (Figure 1b). The study area was divided into three primary geomorphic units and 14 secondary geomorphic units (Figure 1c). Each geomorphic unit presents gradient changes from northwest to southeast in the spatial distribution of geographic elements such as precipitation, temperature, elevation, slope, vegetation, and human activity intensity (

Table A1. Basic information of each geomorphic unit.

Primary Geomorphic Unit	Secondary Geomorphic Unit	Area (km2)	NDVI	Elevation (m)	Slope (°)	POP	GDP	PRE	TEM
Inner Mongolia Plateau	Plateau Hills (I1)	7,376	0.55	1,455	5.35	71	144	432	3.06
	High Plain (I2)	5,722	0.51	1,372	4.05	77	156	401	4.03
North China Mountains	Mountain-Yinshan Eastern Section (II11)	24,192	0.67	1,277	12.48	91	279	494	4.63
	Mountain-Yanshan Section (II12)	29,313	0.79	830	15.09	119	431	579	7.01
	Mountain-Taihang Section (II13)	20,638	0.76	945	18.03	204	727	549	8.42
	Hilly Yanshan Section (II21)	16,687	0.77	262	11.01	329	2,045	616	10.37
	Hilly Taihang Section (II22)	15,561	0.70	238	8.22	481	1,666	499	12.46
	Mountain Basin (II3)	6,930	0.61	515	2.77	366	2,255	515	9.80
North China Plain	Alluvial-Proluvial Fan (III1)	42,450	0.74	30	1.99	1,200	7,696	533	13.09
	Flood Plain (III2)	16,643	0.74	13	2.19	615	1,682	528	13.63
	Yellow River Floodplain (III3)	5,013	0.76	24	2.08	695	1,595	553	13.93
	Depression (III4)	8,544	0.73	11	1.58	695	5,371	538	13.32
	Alluvial and Coast Plain (III5)	10,955	0.68	5	1.51	1073	6,756	552	12.93
	Marine Plain (III6)	4,627	0.47	6	1.27	799	17,389	559	12.74

Notes: NDVI-average NDVI, POP-population density in 2015 (people/km²), GDP-gross domestic product density in 2015 (Yuan/km²), PRE_MN-annual precipitation (mm), and TEM_MN-annual average temperature (℃).

). Due to the policy of national economic development and environmental protection, the ecological environment of BTH has gradually improved [30].

2.2. Data Sources and Preprocessing

Monthly MODND1M China 500 m NDVI products were used as the most important remote sensing data and were downloaded from the Geospatial Data Cloud (http://www.gscloud.cn) and the International Scientific and Technical Data Mirror Site, Computer Network Information Center, Chinese Academy of Sciences. This data set was obtained from MODND1D data using the maximum synthesis method. The MODND1D data were processed by MOD09GA through splicing, cutting, projection conversion, unit conversion, and other processing. The MODND1M data have good reproducibility for MODIS NDVI 16-day synthetic data with 250 m spatial resolution [31]. The Savitzky-Golay filtering method was applied to control the impact of noise in nonphysical data [32]. April to October was identified as the growing season for the major crops [33] for the annual scale NDVI data synthesis. The maximum synthesis method was used to synthesize the NDVI data during the growing season.

To study the spatial differentiation characteristics of the NDVI and analyze the relationship between NDVI trends and different geographical elements, the index selection and processing methods of other researchers were referenced [34,35,36], and the geographic factors were grouped into seven categories and 14 variables (

Table 1. Factors used as explanatory variables for NDVI trends.

Variable Class	Variable Name	Definition and Units	Data Source	Spatial Resolution
Climatic and groundwater	PRE_MN*1	Annual mean precipitation 1980–2015 (mm/yr)	Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences a	1 km
	TEM_MN*1	Annual mean temperature 1980–2015 (℃/yr)		1 km
	PRE_ BT*1	t-test grading of precipitation trends (OLS) during the growing season 2000–2015	China Meteorological Data Network b	1 km
	TEM_ BT*1	t-test grading of temperature trends (OLS) during the growing season 2000–2015		1 km
	GRACE*1	Liquid water equivalent thickness trend (OLS) from GRACE 2003–2015 (cm)	JPL/GRACE-TELLUS c	1°
Geomorphic	ELEVATION*1	Elevation represents macroscopic landform (m)	Geospatial data cloud d	30 m
	SLOPE*1	Slope represents microtopography (°)		30 m
Socioeconomic activities	NLIGHT*2	Nighttime light intensity trend (OLS) from DMSP/OLS 2000–2013	NOAA’s National Geophysical Data Center e	1 km
Accessibility	DIST_RIV*1	Euclidean distance from river (m)	National Catalogue Service for Geographic Information f	vector
	DIST_ROAD*2	Euclidean distance from main road (m)
Agricultural activities	FERT*2	t-test grading of fertilizer uses trend (OLS) 2000–2015	Beijing Municipal Bureau of Statistics g	vector
			Tianjin Bureau of Statistics h
	IRRI*2	t-test grading of effective irrigated area trend (OLS) 2000–2015	Hebei Provincial Bureau of Statistics i
Land cover change	LAND*2	Land cover type change between 2000 and 2015 (6 classes)	Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences j	30 m
Initial NDVI	2000NDVI*3	NDVI in the growing season of 2000	Geospatial data cloud d	500 m

*¹ represents natural variables. *² represents anthropogenic variables. *³ represents a variable with both natural and anthropogenic attributes. ^a http://www.resdc.cn. ^b http://data.cma.cn/. ^c http://upwell.pfeg.noaa.gov. ^d http://www.gscloud.cn. ^e https://www.ngdc.noaa.gov. ^f http://www.webmap.cn. ^g http://www.bjstats.gov.cn. ^h http://stats.tj.gov.cn. ⁱ http://www.hetj.gov.cn. ^j http://www.rcees.ac.cn.

).

The annual precipitation data and annual average temperature data were collected from the Data Center for Resources and Environmental Sciences, the Chinese Academy of Sciences (RESDC). The data set with a 1-km spatial resolution was based on daily observation data from more than 2400 meteorological stations nationwide and was generated by collation, calculation, and spatial interpolation. The monthly temperature and precipitation point data were downloaded from the China Meteorological Data Network. Meteorological observation stations in the BTH and surrounding areas were selected. ANUSPLIN software [37] (Fenner School of Environment and Society at the Australian National University, Linnaeus Way, Acton ACT, Australia) was used to generate monthly raster data with a 1-km spatial resolution, and precipitation and temperature data were synthesized for the growing season.

The composition, distribution, and growth of vegetation are closely related to groundwater, especially shallow groundwater [38]. The Gravity Recovery And Climate Experiment (GRACE) satellite contributed to the information on the distribution of terrestrial water reserve changes by measuring the gravity fluctuations caused by mass anomalies in different regions of the world. This provided a new method for indirectly monitoring NDVI trends [39]. The GFZ GRACE Products (version RL05) (https://grace.jpl.nasa.gov/) from April 2002 to July 2016 were collected, and monthly data for 2003 to 2015 were selected to calculate the trend of the groundwater changes.

As a major indicator of temperature, humidity, and illuminance, ELEVATION was used to characterize the geomorphic gradient effect at the macroscopic scale [40]. SLOPE describes the degree of microscopic changes in surface morphology [25], which was used to reflect the stability of surface material and movement and retention of water.

The Defense Meteorological Satellite Program Operational Linescan System (DMSP/OLS) nighttime imagery data derived from the National Oceanic and Atmospheric Administration (NOAA) and the National Geophysical Data Center (NGDC) revealed information closely related to the distribution of factors such as population and urban area. When using DMSP/OLS light data, socioeconomic factors need not be considered separately [41]. Stable lighting data products were synthesized from 1992 to 2013 global nighttime light time series data. The NLIGHT from 2000 to 2013 was calculated and used as a substitute for the change in human activity intensity over the same period.

Water systems reflect the spatial distribution of surface water and groundwater and control the hydrological conditions of vegetation distribution. Roads indicate human disturbance to natural vegetation and the intensity of the control over artificial vegetation, such as croplands [42]. Data on agricultural facilities, water systems, and roads (including national roads, provincial highways, highways, railways, and urban roads) were collected from the National Catalogue Service for Geographic Information (http://www.webmap.cn). The Euclidean distance calculation method was used to generate DIST_RIV and DIST_ROAD.

Irrigation and fertilization are important for artificially planted crops, and both factors have an impact on crop growth [43]. Data on effective irrigated area and fertilizer use (quantity) from each county except the built-up area in Beijing from 2001 to 2016 were collected from the Economic Yearbook. The annual data were assigned to the administrative unit of each county or district and then converted from vector data into raster data.

Land cover types in BTH include woodland/shrubland, grassland, wetland, cropland, urban, and built-up areas, etc. The dataset was provided by the Research Center for Eco-Environmental Sciences of the Chinese Academy of Sciences (RCEES) (http://www.rcees.ac.cn). Changes in land cover type and the corresponding NDVI trends from 2000 to 2015 were detected, and a numerical reduction in the NDVI was defined as degraded, while an increase was defined as improved. Land cover change was assigned values from −5 to +5 at intervals of 1 depending on the intensity of the change.

The initial NDVI has important effects on the rate, trend, and persistence of NDVI changes, while the factors influencing the changes in the NDVI also have diversity and phase characteristics. The NDVI during the growing season of 2000 was taken as one of the factors influencing the NDVI trends.

2.3. NDVI Trend Detection

The ordinary least squares (OLS) method was used to detect the NDVI trend (θ_slope) [44] from 2000 to 2015, where θ_{slop e} > 0 represents NDVI improvement, θ_slope < 0 represents NDVI degradation, and θ_slope = 0 represents almost no change in the NDVI. To evaluate the consistency of the NDVI trend, a two-sided Student’s t-test for the θ_slope was calculated [45]. Two confidence levels (p < 0.05 and p < 0.01) were identified such as the null hypothesis H0: (θ_slope) = 0. When the confidence level was statistically significant, two types of significance were defined: significant (0.01 < p < 0.05) and extremely significant (p < 0.01).

2.4. Analysis of NDVI Trends Based on GWR

The geographic weighted regression (GWR) model was used to analyze the influence mechanisms of 14 explanatory variables on the NDVI trends. Brunsdon et al. [46] provided a detailed introduction to the principles of GWR, which is expressed in Equation (1).

$$y_i={\beta }_{i0}\left(u_i,v_i\right)+{\displaystyle \sum }_{k=1}^n{\beta }_{ik}\left(u_i,v_i\right)x_{ik}+{\epsilon }_i$$

(1)

where n is the size of the sample, $\left(u_i,v_i\right)$ is the spatial position of the ith sample point, ${\beta }_{i0}\left(u_i,v_i\right)$ and ${\beta }_{ik}\left(u_i,v_i\right)x_{ik}$ are the constant estimates of the ith sample point and the parameter estimate, and ${\epsilon }_i$ (i = 1, 2, 3, …, n) is the random error of the ith sample point. The influence of the observations around the spatial position i on the parameter estimation of the i-point decreases with the increase of the distance, which is represented by the distance attenuation effect. The distance-weighted ordinary least squares method can be used to estimate the parameters below.

$$\widehat{\beta }\left(u_i,v_i\right)={\left(X^{T}W\left(u_i,v_i\right)X\right)}^{-1}{X}^{T}W\left(u_i,v_i\right)Y$$

(2)

where $\widehat{\beta }$ is the estimated value of the parameter β, which is a matrix composed of the independent variable observations. $Y$ is the variable composed of the independent variable observations and W is the spatial weight matrix, which is generally determined by using the Gaussian distance weight method.

$$W\left(u_i,v_i\right)=e^{-\frac{1}{2}}{\left(\frac{d_{ij}}{b}\right)}^2$$

(3)

where $W\left(u_i,v_i\right)$ is the spatial weight of the observation point j for the ith sample point, b is the baseband width of the kernel function, and $d_{ij}$ is the Euclidean distance from the regression point i to the observation point j. When there is no spatial difference between the dependent variable and the independent variable, the data has spatial stability. Brunsdon et al. [47] introduced a stationarity index to evaluate the data stationarity in a regression model.

$$SI=\frac{\beta GWR\_iqr}{2\times GLM\_se}$$

(4)

where $SI$ is the stationarity index, $\beta GWR\_iqr$ is the interquartile range of the standard error of the $x_n$ coefficient, and $GLM\_se$ is the standard error based on the global ordinary least square regression. If $SI\le 1$, the relationship between the dependent variable and the independent variable is stationary, and the smaller the value, the higher the stationarity. The significance test of the model and the spatial autocorrelation of the model error term is performed on the GWR model. AIC [48] was commonly used to test the significance of different models.

$$AIC=2nln\left(\widehat{\sigma }\right)+nln\left(2\pi \right)+n\left[\frac{n+tr\left(S\right)}{n-2-tr\left(S\right)}\right]$$

(5)

where n is the size of the sample, $\widehat{\sigma }$ is the standard deviation of the error term estimate, $tr\left(S\right)$ is the trace of the projection matrix S of the GWR model, which is a function of the bandwidth, and the model with a smaller AIC value is better [49].

In the present study, the GWR model was implemented in the GWR4 software [35,50], which was developed at NCG (National Center for Geocomputation, National University of Ireland Maynooth) and Department of Geography, Ritsumeikan University, Japan (http://gwr.nuim.ie). The GWR4 software was used to complete the regression analysis of the NDVI trends and influencing factors. All the data involved in the calculation were resampled to raster data with a spatial resolution of 1 km, and a point file covering the study area was constructed with each data point controlling nine grids (3 km × 3 km), and the mean value of the nine grids was assigned to the center point. The NDVI change trend was defined as the independent variable, and it was divided into six scenarios: extremely significant decrease areas (abbreviated as NDVI-ESD), significant decrease areas (NDVI-SD), nonsignificant decrease areas (NDVI-ND), nonsignificant increase areas (NDVI-NI), significant increase areas (NDVI-SI), and extremely significant increase areas (NDVI-ESI). These six scenarios were assigned values of −3, −2, −1, 1, 2, and 3. Fourteen factors were defined as the dependent variables, in which the change trends of PRE_BT, TEM_BT, FERT, and IRRI exhibited significant differentiation patterns at the regional scale, so the NDVI trend classifications at the six levels were assigned by these factors. The change trends of GRACE and NLIGHT did not exhibit significant differentiation patterns in the BTH area, and their change trends were still used as dependent variables. In the GWR4 software, the local regression method was chosen, and the accuracy was controlled by optimizing the calculation bandwidth. The GWR-estimated coefficients of the 14 independent variables were used to assess the magnitudes of their influence on the NDVI trends.

2.5. Calculation of the Ratio of Human to Natural Factors

To distinguish the combined effects of different types of explanatory variables on the NDVI trends, the ratio of human to natural factors ($Rati{o}_{h/n}$) was constructed.

$$Ratio{n}_{h/n}=\frac{sum\_human}{sum\_natural}=\frac{{\Sigma }_{i=1}^n\left|est\_huma{n}_i\right|}{{\Sigma }_{i=1}^m\left|est\_natura{l}_i\right|}$$

(6)

where n is the number of natural variables, $sum\_human$ is the sum of the absolute values of all estimated human variable coefficients, $est\_natura{l}_i$ is the estimated coefficient of the ith natural variable, m is the number of human variables, $sum\_natural$ is the sum of the absolute values of all estimated natural variable coefficients, and $est\_huma{n}_i$ is the estimated coefficient of the ith human variable. The comprehensive result of both natural and human effects was 2000NDVI. Thus, it was not considered in this $Rati{o}_{h/n}$ calculation.

3. Results

3.1. Spatial Characteristics of NDVI Trends

The average NDVI from 2000 to 2015 generally increased with a NDVI growth rate (abbreviated as NDVIGR) of approximately 0.005/yr in the BTH area. The area with high NDVI values remained unchanged or slightly increased, and the areas with relatively low and median NDVI values showed a constant or decreasing trend. The fastest NDVIGR in the primary geomorphic units occurred in the North China Mountains (0.006/yr), the slowest NDVIGR occurred in the North Plain (0.0044/yr), and the NDVIGR was 0.005/yr in the Inner Mongolia Plateau. In the edge and transition zones of geomorphic units, the improvement in the NDVI was relatively lower than that in the interior of geomorphic units. For example, the NDVIGR of the Alluvial-Proluvial Fan was only 0.0025/yr. Among the six types of land cover, the NDVIGR of grassland was 0.0072/yr, that of woodland/shrubland was 0.0057/yr, and that of cropland was 0.0053/yr. From 2000 to 2015, the average NDVI of BTH was positively correlated with PRE_MN (coefficient = 0.54, p-value = 0.03) and not significantly negatively correlated with TEM_MN (coefficient = −0.30, p-value = 0.25). NDVI showed a significant reduction when precipitation decreased and temperature increased, and this impact lasted until the following year. This scenario occurred in 2002, 2007, 2009, and 2014 (Figure 2).

The NDVI trends were divided into six scenarios (Figure 3), as shown in Section 2.4. The degraded vegetation area accounted for 8.96% of BTH, of which NDVI-ESD accounted for 0.02%, NDVI-SD accounted for 0.75%, and NDVI-ND accounted for 8.19%. The improvement area accounted for 91.04%, of which NDVI-NI accounted for 40.81%, NDVI-SI accounted for 45.94%, and NDVI-ESI accounted for 4.29%. The total area percentages of NDVI-SI and NDVI-ESI in various land cover types were woodland/shrubland (61.78%), grassland (59.16%), cropland (46.20%), other (39.30%), wetland (34.42%), and urban and built-up areas (31.31%). Meanwhile, the geomorphic gradient exhibited a significant differentiation effect on vegetation changes. NDVI degradation mainly occurred on urban and built-up areas and cropland, and high degradation rates in BTH occurred in the Marine Plain (2.91%), Alluvial-Proluvial Fan (2.12%), and Mountain Basin (2.10%) areas. The highest rates of nonsignificant changes in the NDVI occurred in the Alluvial-Proluvial Fan (73.30%), Hilly Taihang Section (71.98%), and High Plain (70.99%) areas, where the main land types were cropland and grassland. In general, the NDVI improvements at the edge and transition zones of geomorphic units were small, and the greatest improvements occurred inside geomorphic units, in which the Alluvial-Proluvial Fan with the highest population density and high intensity of human activity experienced the most complex changes in the NDVI.

3.2. Multivariate Regression of NDVI Trends

The results of GWR showed that each independent variable (Figure A1) had a clearly different interpretation effect (Figure 4). GWR estimated coefficients (abbreviated as GWRec) indicated the correlation of the independent variables with the NDVI changes. The larger the absolute value of GWRec was, the greater the influence was. In general, the initial NDVI had the greatest impact on the NDVI trends, which was followed by socioeconomic activities, geomorphic factors, climatic factors, groundwater, and agricultural activities, land cover changes, and accessibility factors. All independent variables in grassland, woodland/shrubland, and cropland were highly sensitive to NDVI changes.

To further characterize the influences of natural, human, and their integrated effects on NDVI trends, we calculated three combined parameters of $su{m}_{natural}$, $su{m}_{human}$, and $Rati{o}_{h/n}$. The $su{m}_{natural}$ presented an overall east-west differentiation pattern and further showed a north-south difference (Figure 4o), with two high-value areas in Beijing-Zhangjiakou and Baoding-Cangzhou. The largest $su{m}_{natural}$ occurred in Inner Mongolia, and the value in the North China Plain was slightly larger than that in the North China Mountains. The $su{m}_{natural}$ values of each land cover type are ranked as follows: grassland (1.76), other (1.55), urban and built-up areas (1.51), cropland (1.50), woodland/shrubland (1.45), and wetland (1.41).

As shown in Figure 5a in the Inner Mongolia Plateau, the vegetation of all land cover types except for wetlands generally changed from degradation to vegetation improvement as the natural factors weakened. The vegetation in grasslands, woodlands/shrublands, and other land cover types had the largest response to $su{m}_{natural}$, while the vegetation response was smallest in croplands and urban and built-up areas with the same changes. In the North China Mountains, except for woodland/shrubland and cropland, with the increase in $su{m}_{natural}$, the NDVI of all other land cover types were more inclined to improve. The low $su{m}_{natural}$ in the woodland/shrubland was conducive to the improvement of vegetation. The NDVI improvements in croplands had high requirements for natural conditions even though the $su{m}_{natural}$ was significantly weaker than that in the significant improvement area. In the North China Plain, from the NDVI-ESD to NDVI-ESI, the $su{m}_{natural}$ of cropland and wetland increased slightly but decreased slightly after the NDVI-SI. The $su{m}_{natural}$ of urban and built-up areas showed a continuous state of small increases. The woodland/shrubland and other land cover types showed a monotonous decreasing trend. The $su{m}_{natural}$ of grassland first decreased and then increased. The NDVI-ESI had the highest value.

The $su{m}_{human}$ showed an overall north-south differentiation pattern, followed by an east-west difference (Figure 4p). The high-value areas were mainly distributed in Tianjin and Shijiazhuang in the North China Plain, and the low-value areas were mainly located in the northern and northwestern areas of Zhangjiakou and Chengde. The $su{m}_{human}$ values of the different land cover types were ranked as follows: urban and built-up areas and wetland (0.49), cropland (0.47), woodland/shrubland (0.43), other (0.42), and grassland (0.39). As shown in Figure 5b, from the Inner Mongolia Plateau to the North China Plain, the intensity of human activity increased to approximately 1.5 times that of the plateau. Compared with $su{m}_{natural}$, the $su{m}_{human}$ values of the six land cover types generally exhibited the same trends in the three primary geomorphic units, and the difference in $su{m}_{human}$ was small. The vegetation easily improved with the increase in $su{m}_{human}$ in the Inner Mongolia. However, the regions with less human activity were more prone to vegetation improvements in the North China Plain. In the North China Mountains, the vegetation remained unchanged when the intensity of human activity was high. The vegetation degraded with limited or low-intensity human activity, and the human activity was the weakest in the vegetation improvement area.

The $Rati{o}_{h/n}$ value reflects the combined effects of human factors and natural factors and was generally characterized by an east-west differentiation pattern (Figure 4q), which has two high-value areas. The first high-value area occurred in the curved zone of Qinhuangdao-Tianjin-Langfang-Baoding, which forms a diamond-shaped enclosure with the North China Plain-North China Mountains boundary. The second high-value area was roughly parallel to the boundary with the North China Plain-North China Mountains and was distributed in the central and western regions of Shijiazhuang-Xingtai-Handan. The low-value areas were mainly distributed in Zhangjiakou, Beijing, and Northwestern Chengde. The average of $Rati{o}_{h/n}$ with a distinct geomorphological gradient was 0.16 in the Inner Mongolia Plateau, 0.33 in the North China Mountains, and 0.37 in the North China Plain. The $Rati{o}_{h/n}$ values of each land cover type were ranked as follows: wetland (0.38), urban and built-up areas (0.36), cropland (0.34), woodland/shrubland (0.32), other (0.30), and grassland (0.26).

As shown in Figure 5c, a relatively consistent upward trend of six land cover types occurred in the Mongolia Plateau from NDVI-SD to NDVI-ESI, where the grassland was most sensitive to the changes in $Rati{o}_{h/n}$. In the North China Mountains, the $Rati{o}_{h/n}$ of each land cover type showed a decreasing trend from NDVI-ESD to NDVI-ESI even though the improvement of the woodland/shrubland revealed a high dependence of $Rati{o}_{h/n}$. In the North China Plain, the $Rati{o}_{h/n}$ changes presented the most complex trend between the six land cover types. As the $Rati{o}_{h/n}$ decreased, the vegetation in cropland, grassland, and urban and built-up areas easily improved. Vegetation in woodland/shrubland and other land cover types tended to improve when the $Rati{o}_{h/n}$ increased. Degradation and invariance in wetland vegetation occurred at high $Rati{o}_{h/n}$ values, whereas vegetation increased as $Rati{o}_{h/n}$ increased with reasonably controllable human activities.

3.3. Regression Differences in Different Geomorphic Units

The Inner Mongolia Plateau exhibited the largest $su{m}_{natural}$ and the lowest $Rati{o}_{h/n}$ in BTH, and the NDVI was easily improved when the $Rati{o}_{h/n}$ was high. Due to the conversion of cropland and grassland to urban and built-up areas, NDVI-SD was scattered in the transition zone between the High Plain and the Plateau Hills. NDVI-ND was mainly distributed in the north-central part of this area where cropland was transformed to urban and built-up areas. NDVI-NI largely occurred in cropland and grassland areas and was especially high in the High Plain areas. It was mainly affected by elevation, temperature, precipitation, and agricultural activities. NDVI-SI was mainly distributed in the southeastern part of the Plateau Hills, and it was influenced by geomorphic effect, precipitation, human activities, and warming. The NDVI-ESI trends were generally consistent with that of the NDVI-SI and with the improvement of agricultural technology, and it was mainly distributed in the transition zone adjacent to the North China Mountains.

In the North China Mountains, vegetation changes were mainly driven by natural factors. NDVI-ESD and NDVI-SD were affected by urban-rural expansion and facility construction, and were mainly distributed in the Hilly Yanshan Section, mountain basins, and river valleys. NDVI-ND was mainly distributed in the transition zone of geomorphic units, and affected by urban expansion, infrastructure construction, and precipitation. NDVI-NI was highly correlated with annual mean precipitation and interannual temperature changes and mainly distributed in the Hilly Taihang Section. NDVI-SI was mainly affected by the initial NDVI, elevation, human activities (e.g., returning farmland to forests, small watershed management, ecological economy development), temperature, and precipitation. It was widely distributed in mountains and hills, which was the largest contiguous area of vegetation improvements in BTH and was also the main area of the Three-North Protection Forest Project and the Beijing-Tianjin Sand Source Control Program in BTH. NDVI-ESI occurred in the key area of the Three-North Shelterbelt and was mainly distributed in Chengde and Zhangjiakou. Compared with NDVI-SI, NDVI-ESI had a better combination of precipitation, groundwater volume change, and temperature. Human activities had a greater impact.

In the North China Plain, which had the strongest human activity and the largest $Rati{o}_{h/n}$, the high $Rati{o}_{h/n}$ values corresponded to nonsignificant changes and degradation in vegetation, and the low values indicated vegetation improvement. The vegetation-degraded areas were mainly distributed in the marginal expansion areas of megacities. The NDVI-ND was mainly distributed in the east of the Alluvial-Proluvial Fan and the Marine Plain, where the initial NDVI, agricultural activities in cropland, construction activities in urban-rural areas, and precipitation mainly led to nonsignificant degradation in vegetation. Compared with the NDVI-ND, groundwater, slope, and temperature had prominent effects on the NDVI-NI, which was distributed over large areas of cropland, urban and built-up areas, and wetland. The NDVI-SI was mainly distributed in the eastern part of the Alluvial-Proluvial Fan, where the cropland was mainly traditional low-yield and medium-yield fields and influenced by the initial NDVI, socioeconomic activities, precipitation, temperature, groundwater, and fertilization. NDVI-ESI was mainly distributed in Eastern Hengshui, Southeastern Cangzhou, and Northern Handan. The roles of agricultural production and management technology in cropland and greening projects in urban and built-up areas were more prominent. The traditional low-medium-yield fields were transformed into medium-high-yield fields, which promoted the extremely significant increase in vegetation. In general, vegetation improvements exhibited the characteristics of regional concentration and convergence of driving factors and $Rati{o}_{h/n}$, while vegetation degradation had many influencing factors and strong heterogeneity in the North China Plain.

4. Discussion

4.1. Impact of Climate on NDVI Trends

The NDVIGR of this study is 0.005/year, which is higher than the NOAA-GIMMS-NDVI in the area [51,52] and is basically consistent with the SPOT-VGT-NDVI [28,53] and previous findings using the MODIS NDVI [54,55]. The fluctuations in the NDVI throughout the region are closely related to meteorological events such as droughts. The sharp decline in the NDVI in 2009 was attributed to the impacts of the 2008 winter and 2009 drought [56]. The continuous decline in the NDVI from 2014 to 2015 was affected by drought events in the Northern Hemisphere in 2014 and 2015 [57,58].

Since climate factors have long-term and cumulative effects and are highly integrated with human activities, precipitation and temperature do not always have a decisive influence on the interpretation of the causes of NDVI changes. From 2000 to 2015, the sensitivity between NDVI trends and PRE_MN was stronger than that between NDVI trends and TEM_MN [55,59]. TEM_BT was more sensitive than PRE_BT. Except for the Mountain-Yinshan Eastern Section, the rest of the BTH area exhibited a warming pattern [60]. The sensitivity of TEM_MN to NDVI changes was great when TEM_BT was low, and there was a certain reconciliation between them [61]. Except for the Marine Plain, PRE_MN declined in each geomorphic unit. The correlation between the NDVI trends and precipitation and temperature exhibited significant spatial differences. The areas with large geomorphological gradients and climatic gradients, such as the low NDVI value areas in plateaus, plain croplands, and grasslands, were sensitive to climatic factors. The NDVI in mountainous forests was less sensitive to changes in climatic factors, such as the >600 mm curved rainy belts distributed on the windward slopes of the Yanshan and Taihang Mountains [62]. The NDVI in the north and northeast of the relatively humid mountainous area of BTH did not always respond directly and proportionally to precipitation changes. Although a nonlinear change in the soil water supply was caused by the increase in precipitation, the water utilization rate of vegetation changed with the changes in temperature and humidity [63]. In general, the NDVI trends were mainly positively and not significantly correlated with PRE_MN or TEM_MN, which accounts for 94.02% of the total area, while the significant and extremely significant positive correlations accounted for only 5.98% of the area, which were mainly distributed in mountain basins and depressions as well as the border area between the Mountain-Yinshan Eastern Section and the Mountain-Yanshan Section. The variation in effective groundwater depth exhibited coupling and hysteresis with precipitation. In the area with extensive groundwater exploitation in the North China Plain [64], especially in the funnel area, the influence of the change in groundwater depth on the NDVI trends correspondingly weakened [65]. Considering the comprehensive effects of climatic factors on NDVI trends, climatic factors had no significant or negative effects on NDVI, which indicates that climate change had a certain negative effect on vegetation changes [66], while human activities had a great impact.

4.2. Influences of Non-Climatic Factors on NDVI Trends

4.2.1. Geomorphology

Both macroscopic and microscopic geomorphology had important impacts on the distribution of soil and hydrothermal conditions, as well as strong coupling effects including the intensity of human activity, land cover, vegetation distribution patterns, and its changes [67]. Land use regulations under various slopes in the National Soil and Water Conservation Law directly affect the type of vegetation cover. High elevations and steep slopes (>25°) in the Inner Mongolia Plateau and North China Mountains are highly suitable for ecological conservation, such as the “Grain to Green Project” and the “Natural Forest Conservation Program” [2]. However, low altitudes (<20 m) correspond to improved water sources, and steep slopes (>4°) with better drainage conditions in the North China Plain are conducive to cropland and woodland vegetation improvements. The combined effect of elevation and slope on the NDVI trends would be further magnified when the internal geomorphic gradient increases, such as in the six land cover types of the High Plateau, the cropland and urban and built-up areas of the mountain basins in the North China Mountains, and the wetland in the depression of the North China Plain.

4.2.2. Land Use

Land use is one of the most important factors in explaining the regional differences in the NDVI attributes that were analyzed. The transformation of land cover type has a direct impact on NDVI [68] in only the occurrence area, and the influences are limited in the areas at the edge and outside of the boundaries. From 2000 to 2015, land cover type conversion occurred in 7.5% of BTH, with the highest sensitivity to land conversion in the plain wetlands. The degraded NDVI areas were mainly distributed within 200 to 500 m from the urban and built-up areas, and the improvement areas were distributed by 1 km or more from the urban and built-up areas. Within 1 km of the land cover type change area, the NDVI stability or improvement trends were suppressed, and the intensity of this influence was gradually attenuated outside 3 km. The NDVI improvements caused by land cover changes mostly occurred in mountainous and hilly areas, while, in plains and mountain basins, land cover transformation mainly led to nonsignificant changes in vegetation.

4.2.3. Accessibility

Roads and rivers reflect the different needs of people and farming. Vegetation degradation is concentrated in the residential areas near water and roads. The roads used in the present study are mainly high-grade roads, which mainly inhibit vegetation improvement, while low-grade roads closely related to cropland are not all involved in the calculation. The water system in BTH is highly influenced by human activities, such as the construction of reservoirs, water diversion, and river regulation. Thus, the morphological structure of the water system does not fully represent the water volume distribution. Urban and industrial production and agriculture have intensified water stress, especially during the growing season. Surface water irrigation and groundwater irrigation have become the main countermeasures in the North China Plain [69], which resulted in vegetation changes associated with off-site water and groundwater. Therefore, the impact of the water system on local vegetation changes is weakened. The water system distributed in the mountainous area is relatively closely related to the vegetation of the catchment basin, which is associated with terrain, elevation, slope, and aspect [70]. Vegetation near the water system can be easily improved even though the effect of the water system may also be weakened by the terrain effect [71].

4.2.4. Agricultural Production Technology

Cropland accounts for 42% of BTH, which is the largest land use type. The low GWRec associated with the irrigation and fertilization practices in cropland reflect the direct impact of agricultural activities on vegetation changes and reveal the indirect effects of other agricultural activities on vegetation changes, such as soil improvement, seed improvement, new technology use, and planting structure improvement. These activities are used to actively adapt to climate change and regional economic development and have become important factors affecting cropland vegetation [72]. In the medium-yield and low-yield fields, vegetation has been significantly increased through soil improvement, improving irrigation and fertilization efficiency, and using other technologies and management measures. Except for the Inner Mongolia Plateau and the Southeastern part of the North China Plain, the effects of irrigation and fertilization on increasing yields have been reduced, especially in the medium-high-yield fields [73].

4.2.5. Socioeconomic Activities

This variable is negatively correlated with NDVI trends, which suggests that activities such as urban-rural spatial expansion, spillover effects, and impervious surfaces may inhibit vegetation improvements [74]. The statistics also show that the growth rate of NLIGHT in the plain vegetation improvement area is only one-half or one-third of that in the vegetation degradation area, and the rate decreases in mountains and plateaus. However, in cropland, ecological project areas, and large cities, human activities can effectively promote vegetation improvement by addressing climate and geomorphological factors. For example, the government has implemented ecological economic development measures in forest and grassland and improved agricultural infrastructure in cropland. These measures have promoted the adjustment of urban-rural industrial structure, and more optimized human activities have accelerated the improvement of vegetation. The improvement of vegetation in megacities is affected by many effects, such as urban green space improvement, urban heat island effects, increased CO₂ concentrations, and increased active nitrides, all of which are highly correlated with human activities [75]. However, the intensity of human activity may be affected by the saturation effect of DMSP/OLS light data, especially in the bright cores of urban centers [76].

4.2.6. Initial NDVI

The 2000NDVI represents the initial state of vegetation and constrains the derived information, such as the direction, strength, and reliability. The average estimated coefficient of BTH is approximately −0.36, which is roughly similar to geomorphology and human activity factors. The negative correlation between the 2000NDVI and NDVI trends reflects two points of information: first, there may be a saturation effect on NDVI data as well as DMSP/OLS lighting data, especially in the high NDVI value areas, and, second, the extremely low and very high-value regions of vegetation are inversely changed to some extent [77]. For example, low-coverage vegetation was improved, such as various land cover types in the Inner Mongolia Plateau, and high vegetation cover was degraded, such as woodland/shrubland and grassland in the North China Plain and woodland/shrubland in the North China Mountains.

4.3. Implications and Limitations

Following the gradient of climate and geomorphology, people should actively respond to changes in the natural environment as well as regulate and optimize the layout and intensity of human activities. The stability of vegetation in the North China Mountains should be maintained in the future. The intensity of human activities in the Inner Mongolia Plateau should be reasonably controlled and the planting structure and agglomerated woodland/shrubland should be optimized to effectively control the scale of negative effects from urban and built-up areas in the North China Plain. Based on functional ecological zoning, suitable ecological economic developments should be preferred in different functional zoning units and the impact of human activities on the ecological environment should be scientifically guided and controlled.

The spatial resolutions of the 14 interpretation factors used in this study are different and do not match the spatial resolution of NDVI, such as groundwater data, irrigation data, and fertilization data. Meanwhile, because of the saturation effect of NDVI, the maximum synthesis method may reduce data accuracy at the interannual scale, which leads to unsustainable or noncontinuous conditions in NDVI trends, especially those affecting the evaluation of the NDVI change trends. Although the impacts of human activities were characterized by using nighttime light intensity, roads, water systems, irrigation, fertilization, and other indicators, there is still a gap in the expression of diversity, complexity, and dynamics of human activities, especially in cropland areas. Although collinearity is generally evaluated, there may be individual factors in the interpretation of the GWR, which are repeated in the local space and lead to the repeated or amplified interpretation of certain factors.

5. Conclusions

This study shows that the overall vegetation in the BTH region has improved over the last 16 years, and NDVI fluctuations in individual years were closely related to climate events. Vegetation changes showed a northwest-southeast pattern, which was highly consistent with the geomorphic gradient. The main vegetation improvement areas were distributed in the North China Mountains where natural effects were pronounced, and the main vegetation invariance and degradation areas were distributed in the North China Plain where human effects were prevalent. Natural factors were generally more sensitive than human factors to vegetation changes, and the $Rati{o}_{h/n}$ better reflected the combined effects of human and natural factors. The Inner Mongolia Plateau had the largest $su{m}_{natural}$ and the lowest $Rati{o}_{h/n}$ in the BTH, and the vegetation was more inclined to improve in areas with high $Rati{o}_{h/n}$ values. In the North China Mountains, the high $Rati{o}_{h/n}$ corresponded to vegetation degradation and invariance, while the low values implied vegetation restoration due to the implementation of ecological engineering. The vegetation improvements had the characteristics of regional concentration and convergence of driving factors and $Rati{o}_{h/n}$, while vegetation degradation had many influencing factors and strong heterogeneity in the North China Plain. The results of this research can assist with ecological functional zoning and ecological economic development in BTH. Future related research should further refine the impacts of human activities, improve data accuracy, and reduce the saturation of NDVI data.