Spatial–Temporal Vegetation Dynamics and Their Relationships with Climatic, Anthropogenic, and Hydrological Factors in the Amur River Basin

Abstract

Information about the growth, productivity, and distribution of vegetation, which are highly relied on and sensitive to natural and anthropogenic factors, is essential for agricultural production management and eco-environmental sustainability in the Amur River Basin (ARB). In this paper, the spatial–temporal trends of vegetation dynamics were analyzed at the pixel scale in the ARB for the period of 1982–2013 using remotely sensed data of long-term leaf area index (LAI), fractional vegetation cover (FVC), and terrestrial gross primary productivity (GPP). The spatial autocorrelation characteristics of the vegetation indexes were further explored with global and local Moran’s I techniques. The spatial–temporal relationships between vegetation and climatic factors, land use/cover types and hydrological variables in the ARB were determined using a geographical and temporal weighted regression (GTWR) model based on the observed meteorological data, remotely sensed vegetation information, while the simulated hydrological variables were determined with the soil and water assessment tool (SWAT) model. The results suggest that the variation in area-average annual FVC was significant with an increase rate of 0.0004/year, and LAI, FVC, and GPP all exhibited strong spatial heterogeneity trends in the ARB. For LAI and FVC, the most significant changes in local spatial autocorrelation were recognized over the Sanjiang Plain, and the low–low agglomeration in the Sanjiang Plain decreased continuously. The GTWR model results indicate that natural and anthropogenic factors jointly took effect and interacted with each other to affect the vegetated regime of the region. The decrease in the impact of precipitation to vegetation growth over the Songnen Plain was determined as having started around 1991, which was most likely attributed to dramatic changes in water use styles induced by local land use changes, and corresponded to the negative correlation between pasture areas and vegetation indexes during the same period. The analysis results presented in this paper can provide vital information to decision-makers for use in managing vegetation resources.

1. Introduction

As a major component of terrestrial ecosystems, vegetation plays an important role in material cycling and energy flows, and provides irreplaceable service functions that maintain the wellbeing of our planet and all the creatures that inhabit it. These function services include food provision, climate regulation, carbon sequestration, timber production, biodiversity preservation, and soil protection [1,2,3,4,5,6,7,8]. Vegetation growth affects the ecological balance, the terrestrial carbon cycle, water circulation, and other biochemical processes [2,9]. Thus, information about the spatial–temporal dynamics of regional and global vegetation is a fundamental need in order to facilitate the better management of our planet ecosystem, provide effective support for environmental sustainability and ensure the safety of agricultural production; these have become emerging issues in the field of environmental studies [2,5,10].

With the rapid advancement of remote sensing technology, global-scale land surface parameters, derived from long time series of remote sensing observations, have provided a good description of vegetation changes [11,12,13]. Leaf area index (LAI), defined as one-half the total developed area of green leaves per unit horizontal surface area, can serve as a proxy for detecting vegetation greenness and serves as a measure of the amount of vegetation that is vertically distributed [14,15,16,17,18,19]. Fractional vegetation cover (FVC), defined as the projected area of aboveground vegetation per unit ground area, can characterize the horizontal density of land surface vegetation and reflect plant community structure [20,21,22]. In addition, gross primary productivity (GPP) is the photosynthetic carbon assimilation by land plants per unit space and time, and can serve as an indicator of vegetation productivity [23,24,25]. Based on remote sensing products of these parameters at regional and global scales, scholars around the world have used different methods to study spatial–temporal changes in vegetation greenness (related to canopy structure) and vegetation productivity (related to canopy function) [26,27,28], such as trend analysis and spatial autocorrelation analysis.

Vegetation growth, production and distribution are highly reliant on—and sensitive to—natural and anthropogenic factors, such as precipitation, temperature, water resources, farming, and urban expansion. The condition of the global vegetation is constantly changing at various spatial and temporal scales, driven by natural and anthropogenic factors [29,30,31]. The response of vegetation dynamics to differences in environmental factors varies significantly across regions, according to the regional climate conditions, water availability, and land cover [32,33].

As the tenth largest river basin in the world, the natural conditions of the Amur River Basin (ARB), a transboundary basin in eastern Eurasia, vary greatly from region to region. In addition, the impact intensity of human activities on the ground vegetation has varied considerably under the influence of different national policies, lifestyles, and economic development [34,35,36,37,38]. In recent decades, significant climate changes have been observed in the ARB, where an annual mean temperature warming rate of about 0.34 °C/10a for 1975–2004 was reported, which is approximately 2.5 times higher than that for the 1891–2004 period [39]. Regional climatic characteristics and dramatic climate change govern the spatial–temporal changes of vegetation in the basin. In the context of changes in climatic and hydrological factors, the continued intensification of agricultural activities, driven by the large demand for food and economic development, affects the vegetation dynamics in the ARB and threatens the ecosystem [40,41,42].

Most previous studies have used correlation analyses and regression models to reveal the factors driving changes in vegetation dynamics in the whole basin or sub-basins of the ARB. For example, the Pearson correlation analysis, applied to investigate the relationship between growing season and seasonal Normalized Difference Vegetation Index (NDVI) and climatic variations in the ARB, was reported in [37] and it was concluded that growing season NDVI was mainly regulated by precipitation. A linear trend analysis of vegetation, temperature, and precipitation, and a comparison between the correlation of NDVI–temperature and NDVI–precipitation in northeastern China was carried out by Mao, et al. [43]. However, the analysis methods used in most of the previous studies focused solely on a temporal scale or spatial scale, and therefore, they do not provide a comprehensive analysis to deal with variability between vegetation dynamics and environmental factors at both spatial and temporal scales simultaneously [44,45,46]. Additionally, a few efforts have been made to investigate the relationships between hydrological factors and vegetation dynamics in the ARB, with close interactions between hydrological factors and vegetation growth having been found in many studies [47,48,49]. Therefore, exploring the spatial–temporal vegetation dynamics and their relationships with climatic, anthropogenic, and hydrological factors in the ARB is essential for better understanding the mechanisms driving vegetation changes, and can provide more detailed scientific support for environmental sustainability and agricultural production management. The objectives of this study are specified as follows:

• Analyze the spatial trends of three parameters (LAI, FVC and GPP) that can be considered as representative of the growth condition of surface vegetation in the study region using Mann–Kendall and Sen’s slope methods to understand the spatial variability of vegetation growth conditions.
• Explore the spatial autocorrelation characteristics of vegetation indexes in the ARB, based on the results derived from Moran’s I technique on remotely sensed vegetation information, for determining the rapid shift of eco-system changes over previous decades.
• Utilize partial least squares regression (PLSR) and geographical and temporal weighted regression (GTWR) models to further evaluate the relationships between land surface parameters and climatic factors, land use/cover types, and hydrological variables simulated with the soil and water assessment tool (SWAT) model in the ARB to clarify spatial–temporal vegetation dynamics and their relationships with climatic, anthropogenic, and hydrological factors in the Amur River Basin.

2. Study Area

As the 10th largest river basin in the world, the ARB (41°–56° north and 107°–142° east) is situated in northeastern Asia. Its tributary rivers flow through China, Russia, Mongolia, and North Korea [36,38,41], covering a total land area of approximately 2 million square kilometers with elevations ranging from 0 (the mouth of the Amur River into the Tartary Strait) to 2565 m asl (the Khentii Mountains in northeastern Mongolia) [50,51] (See Figure 1). As one of the most important transnational river basins, the ARB exhibits quite different regional geographic and climatic characteristics: the topography is characterized by high mountains over the west, with typical continental climate features and low hills in east margin of the watershed, with an extensive plain in between, where the monsoonal climate is predominant. Population density presents a north–south gradient over the basin, with notable different anthropogenic effects [34,35,38]. Obtaining a better understanding of the vegetation dynamics driven jointly by natural and anthropogenic factors in the ARB, is therefore necessary for environmental sustainability studies. The hydrological variables derived from SWAT model simulations were confined in the sub-basin of the ARB where the Komsomolsk hydrological station gauged, covering 94% area of the whole ARB.

3. Materials and Methods

3.1. Datasets

In this study, Global land surface satellite (GLASS) LAI, FVC and GPP products for the period of 1982–2013 were used to analyze vegetation dynamics. These GLASS products were released using the Hierarchical Data Format (HDF), accessible to the public for free at http://glass-product.bnu.edu.cn. The detailed information of land surface parameter products, representative of the vegetation dynamics used in this study are presented in

Table 1. Land surface parameter products used in this study as representative of vegetation dynamics.

Parameter	Product	Sensor	Spatial Resolution	Temporal Resolution	Data processing	Reference
LAI	GLASS	AVHRR	about 5 km	8-day	(1) Use MVC technique to produce monthly data based on 8-day products (2) Average monthly data to obtain annual data	[19,58]
FVC	GLASS	AVHRR	about 5 km	8-day	(1) Use MVC technique to produce monthly data based on 8-day products (2) Average monthly data to obtain annual data	[59]
GPP	GLASS	AVHRR	about 5 km	8-day	Summing 8-day products to obtain annual data	[60]

NOTE: LAI: Leaf area index; FVC: Fractional vegetation cover; GPP: Gross primary productivity; GLASS: Global land surface satellite; AVHRR: Advanced Very High Resolution Radiometer; MVC: maximum value composite.

. The monthly LAI and FVC values were computed using the maximum value composite (MVC) technique [52], which is widely used in compositing NDVI and other vegetation indexes in order to minimize the interference of atmospheric effects, the scan angle, cloud contamination, and solar zenith angle [29,53,54,55,56,57].

According to characteristics of climate changes during different periods in the ARB that have described in previous studies [61], in addition to the land use/cover transformations reported in the literature [36,40,62], the entire study period was divided into four periods to analyze the spatial–temporal relationships between environmental factors and vegetation dynamics: 1982–1990 (Period 1), 1991–1999 (Period 2), 2000–2006 (Period 3), and 2007–2013 (Period 4). The detailed information of datasets used in this study as representative of environmental factors are presented in

Table 2. Datasets used in this study as representative of environmental factors.

Data Type	Description	Source	Download Site
Climate	Daily temperature and precipitation data from 1982 to 2013 in 193 meteorological stations	China Meteorological Data Network (CMA) and National Oceanic and Atmospheric Administration (NOAA)’s National Centers for Environmental Information (NCEI)	https://data.cma.cn/en https://www.ncdc.noaa.gov/cdo-web/
Land-use	Chinese land-use map (1980, 1995, 2005, and 2010) and global land-use map (1992–1993, 2000, 2005, and 2010) at spatial resolution of 1 km	Chinese land-use maps from Resource and Environment Data Cloud Platform; Moderate Resolution Imaging Spectroradiometer (MODIS) land cover type product from U.S. Geological Survey (USGS); Global Land Cover Characterization (GLCC) from USGS; and Global Land Cover 2000 database (GLC2000) from Joint Research Centre, European Commission	http://www.resdc.cn/Datalist1.aspx?FieldTyepID=1,3 https://lpdaac.usgs.gov/products/mcd12q1v006/ https://www.usgs.gov/centers/eros/science/usgs-eros-archive-land-cover-products-global-land-cover-characterization-glcc http://forobs.jrc.ec.europa.eu/products/glc2000/data_access.php
Hydrology	Surface runoff (SURQ), lateral flow (LATQ), snowmelt (SM), soil water (SW), ground-water flow (GWQ), and evapotranspiration (ET)	Simulated by the soil and water assessment tool (SWAT) model	authors’ previous study [61]

. To investigate the driving factors of climate on vegetation changes, temperature and precipitation data observed from 1982 to 2013 in 193 meteorological stations within the basin were used. The climate data were interpolated at a spatial resolution of 8.8 × 8.8 km by using thin plate smoothing splines and Anusplin software, and the LAI, GPP, and FVC datasets in the same duration were resampled to the same spatial resolution. The combined datasets of Chinese land use maps with global land use maps from the closest dates in the four historical periods were used as inputs for the GTWR model to investigate the anthropogenic driving factors responsible for regional vegetation changes (urbanization and wetland reclamation, etc.). In addition, the hydrological variables simulated by the SWAT model, as presented in our previous study [61], were adopted to analyze the relationship between hydrological variables and vegetation dynamics.

3.2. Methods

The paper presents the results of applying the Mann–Kendall, Sen’s Slope, and Global and Local Moran’s I techniques to analyze vegetation changes and applying PLSR and GTWR models to analyze the relationships between vegetation and climatic, anthropogenic, and hydrological factors based on meteorological data, remotely sensed information and hydrological variables simulated by the SWAT model. Figure 2 shows an overview of the general workflow followed in this study.

3.2.1. Time Trend Analysis

The magnitude of the trends in annual LAI, FVC and GPP were estimated using Sen’s slope method [63]. The median value of the slope series computed from two consecutive points of the series can help reduce the influence of outliers or missing data [63,64,65]. The formula of Sen’s slope method is as follows:

$$Slope=\mathrm{Median}\left[\frac{(x_j-x_i)}{(j-i)}\right], \forall j > i$$

(1)

where $x_i$ and $x_j$ are the values at times i and j ($1\le i<j\le n$), respectively, and a positive (negative) value of Slope indicates an upward (downward) trend.

The Mann–Kendall [66,67] non-parametric statistical test, which is frequently used to examine the significance of trends in time series of data, remining insensitive to outliers [23,68,69], was used to analyze the significance of the annual LAI, FVC, and GPP trends in each grid for the ARB during the 1982–2013 period. At a given significance level a = 0.05, the threshold of the normal distribution is $Z_{1-a/2}$ ($Z_{1-a/2}=Z_{0.975}=1.96$). When $\left|Z\right|\le Z_{1-a/2}$, the null hypothesis can be accepted (the trend is nonsignificant), and when $\left|Z\right|>Z_{1-a/2}$, the trend is significant. The sequential version of the Mann–Kendall test proposed by Sneyers [68] was used to test abrupt changes in data series. If the intersection point of the forward sequence statistic UFk and the backward sequence UBk, which are calculated using the same equation but in reverse data series order, is between the boundary lines of the confidence zone ($\pm 1.96$ at the confidence level of 95%), the critical point of change can be detected [70,71,72,73].

3.2.2. Spatial Autocorrelation Analysis

Spatial autocorrelation is the term used to describe the dependency and covariance of variables within a geographic area, and it reveals the phenomenon of geographic proximity’s potential [74,75]. In this paper, both global and local measures of spatial autocorrelation were calculated using the Python package GeoRasters to evaluate the degree of spatial autocorrelation in vegetation growth patterns across the study region and to indicate the variations in spatial autocorrelation throughout the study period [12,76].

Moran’s I [77] is one of the most commonly used indexes of global spatial autocorrelation. Moran’s I values range from −1 to 1: the index value from 0 to 1 indicates a positive autocorrelation, which implies a clustering state of spatial geographical phenomena, and the index value from −1 to 0 represents a negative autocorrelation, which implies a larger difference between neighbors, [74,78,79,80]. The formula for Moran’s I is given in Equation (2):

$$I=\frac{n{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{w}_{ij}\left(x_i-\overline{x}\right)\left(x_j-\overline{x}\right)}}}{{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{w}_{ij}{\displaystyle \sum _{i=1}^n{\left(x_i-\overline{x}\right)}^2}}}}$$

(2)

where $x_i$ and $x_j$ are the values of variable x (vegetation indexes) at location i and location j, respectively. $\overline{x}$ is the average value of x, n is the total number of locations, and $w_{ij}$ represents the spatial weight [4,44,81,82].

A local spatial autocorrelation analysis was performed using the local indicators of spatial association (LISA) [83,84] to indicate the spatial heterogeneity of the geographic phenomena at the local level across the study area, based on the following equation [74,83,85]:

$$I_{\mathrm{i}}=\frac{(x_i-\overline{x})}{S^2}{\displaystyle \sum _j{w}_{ij}(x_j-\overline{x})}$$

(3)

where $S^2$ represents the variance of the observed values. In this study, the LISA map shows the spatial distribution (clustered/random/dispersed) of each vegetation index based on 9999 permutations at the significance level of p < 0.05. The spatial cluster/outlier characteristics of regions can be divided into four categories: low–low clusters (LL), high–high clusters (HH), low–high outliers (LH), and high–low outliers (HL). Positive and negative LISA values represent the aggregation of similar observations—“spatial cluster”, and different observations—“spatial outlier”, respectively [79,80,86,87,88,89].

3.2.3. GTWR

GTWR [45] was used in this study to reveal the spatial–temporal heterogeneity of the relationships between vegetation growth (LAI, FVC and GPP) and environmental variables (climatic, anthropogenic and hydrological factors) in the ARB at the sub-basin level [90]. Unlike the traditional geographically weighted regression model which considers only the spatial dimension, the GTWR model can well reflect evolution in the spatial–temporal relationships between vegetation dynamics and environmental variables and can simultaneously explore both spatial and temporal non-stationarity in vegetation changes [91,92]. The GTWR model can be expressed as follows:

$$Y_i={\beta }_0(u_i,v_i,t_i)+{\displaystyle \sum _k{\beta }_k(u_i,v_i,t_i)}{X}_{ik}+{\epsilon }_i$$

(4)

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

(5)

where $(u_i,v_i,t_i)$ denotes the coordinates of location i in space and time; $Y_i$ and $X_{ik}$ are the dependent variable and the kth independent variable, respectively; ${\beta }_0(u_i,v_i,t_i)$ represents the intercept value, ${\beta }_k(u_i,v_i,t_i)$ is a set of parameter values at space–time location i; $\widehat{\beta }(u_i,v_i,t_i)$ is the estimation of the parameters ${\beta }_k(u_i,v_i,t_i)$; $W(u_i,v_i,t_i)$ denotes the spatial–temporal weight matrix, and ${\epsilon }_i$ is the residual error.

4. Results

4.1. Spatial–Temporal Variations in Vegetation in the ARB

To detect the spatial–temporal variations of vegetation in the ARB basin, the Mann–Kendall test and Sen’s slope were employed to examine the trends and abrupt changes in the LAI, FVC and GPP from 1982 to 2013 at the confidence level of 95%. Figure 3a–c presents the Mann–Kendall Z value and Sen’s slope of the area-averaged annual vegetation indexes across the study region. During the study period, the annual FVC in the study region varied from a minimum of 0.327 in 1983 to a maximum of 0.352 in 2013 with a significant increase rate of about 0.0004/year (Z = 3.57, p < 0.05). While, the changes in area-averaged annual LAI in the same duration were not significant (Z = −1.67, p > 0.05), these values varied from a minimum of 1.369 in 2003 to a maximum of 1.563 in 1999, with a slight decreasing trend of −0.0013/year. Overall, an increasing trend at 0.19 gC m⁻² yr⁻² (Z = 0.21, p > 0.05) was found in annual GPP, however, large interannual fluctuations with a maximum value of 900.37 gC m⁻² yr⁻¹ (1988) and a minimum value of 784.05 gC m⁻² yr⁻¹ (2003) were observed. The period averages of LAI, FVC, and GPP for the four periods in Figure 3a–c indicated that the period averages of FVC gradually increased during the four periods, while the period averages of LAI and GPP changed most dramatically between Period 2 and Period 3 (before and after 2000).

The sequential Mann–Kendall tests for the annual LAI with a forward-trend UFk and backward-trend UBk, as illustrated in Figure 3d–f, indicated an obvious abrupt change in LAI in 2001. In addition, a significant abrupt change in GPP occurred in the study area around 2001. Several abrupt changes in FVC occurred from 1999 to 2008, but most of them were insignificant. These results suggested that all three types of indexes describing different vegetation statuses had undergone a rapid change around the year 2001.

Figure 4a–i shows the spatial patterns of pixel-by-pixel Sen’s slope and the Mann–Kendall Z value estimated over the ARB, as well as the significance of the trends computed at a 95% significance level for the time series LAI, FVC and GPP images. Apparent spatial heterogeneity for the changes of LAI from 1982 to 2013 over the ARB, as exhibited in Figure 4a–c, was carefully investigated. The results indicated that about 37.3% of the basin area showed an increasing tendency in terms of LAI, in which approximately 7.7% tended to increase significantly (p < 0.05), where the most obvious increasing trend occurred over the Songnen Plain (an important food and oil production national base of China) in the southern part of the ARB. On the contrary, a decreasing trend in LAI was significant over about 18.9% of the ARB area, where the most significant reductions in LAI were mainly distributed in Mongolia and the Inner Mongolia Autonomous Region of China in the southwestern ARB area. Consistent with the area-averaged LAI trend, the LAI changes at the regional scale were insignificant over a large portion of the study area.

Over past 32 years in the ARB, 84.3% of the total basin area has exhibited an increasing trend in FVC, and pixels with a significant increasing trend accounted for about 62.5% of this value (p < 0.05), mainly being concentrated in the Chinese and Russian parts of the ARB. Only 2.0% of the total area showed a significantly decreasing trend, and pixels with a decreasing trend were mainly distributed in the southwestern margin of the ARB. For GPP, approximately 45.1% of the area in the ARB presented a decreasing trend, in which about 10.5% of the area was decreasing significantly (p < 0.05). Furthermore, areas with an increasing tendency accounted for 54.9% of the area, among which, 16.3% showed a significant increasing trend over the ARB from 1982 to 2013 (slope > 0, p < 0.05).

4.2. Spatial Autocorrelation Analysis of Vegetation Dynamics in the ARB

The global Moran’s I was investigated for the four period-averaged LAI, FVC, and GPP during 1982–2013 over the ARB. The global Moran’s I values of LAI for the four periods were 0.898, 0.897, 0.913 and 0.917 with p-values all being less than 0.01, respectively, which indicated a significant positive spatial autocorrelation and persistent high–high or low–low aggregations of the LAI. Similarly, the global Moran’s I values of FVC and GPP for the four stages were also all above 0.8, with p-values all less than 0.01, indicating positive spatial autocorrelation and persistent high–high or low–low aggregations. However, the global Moran’s I values do not indicate the possible spatial changes of spatial autocorrelation of the studied indexes across the basin.

Figure 5 presents the clustering maps of the local Moran’s I of LAI, FVC and GPP in the ARB for the four studied periods (Figure 5a–d show Periods 1–4, respectively). The local Moran’s I for LAI in the four periods ranged from −1.83 to −10.41, −1.82 to 9.47, −1.57 to 5.55 and −0.90 to 5.93, respectively. With respect to the local Moran’s I of the LAI dataset for Periods 1–4, the areas of low–low agglomeration were always concentrated in the Mongolian region and the Songnen Plain and the Inner Mongolian Autonomous Region of China over four different phases. The clustering maps of the local Moran’s I values further revealed that the local spatial autocorrelation of the Sanjiang Plain demonstrated the most significant change and the low–low agglomeration in the Sanjiang Plain decreased persistently. Among these changes, the most significant changes were found between Period 2 and Period 3: i.e., the autocorrelation changed the most around 2000 over the Sanjiang Plain, while the low–low agglomeration in the Sanjiang Plain disappeared after 2000.

Figure 5e–h presented the clustering maps of local Moran’s I of FVC in the ARB for four studied periods. The local Moran’s I of FVC for the four periods ranged from −0.67 to 6.45, −0.69 to 6.67, 0.57 to 6.21 and −0.65 to 6.63, respectively. Like the LAI spatial autocorrelation, the local spatial autocorrelations of FVC for the four periods were relatively stable, especially for the high–high agglomeration in the northern part of the basin. The clustering maps for the four periods all displayed strong high–high and low–low agglomerations, whereas the spatial outliers (high–low clustering or low–high clustering) were not conspicuous. A similar decreasing trend in low–low agglomeration areas in the Sanjiang Plain was also found during the study period.

Figure 5i–l presents the clustering maps of local Moran’s I of GPP in the ARB for four studied periods. The local Moran’s I of GPP for the four periods ranged from −1.77 to 8.24, −1.94 to 8.59, −1.70 to 6.27 and −1.90 to 6.93, respectively. For the four periods, the areas characterized with low–low agglomeration were always found in the Mongolia and the Songnen Plain and the Inner Mongolian Autonomous Region of China. Unlike that of the other two indexes, the spatial autocorrelation of GPP on the Sanjiang Plain changed very slightly, with no obvious low–low clustering in this area. The most significant changes in the clustering of GPP were found in Mohe County, the far most northern part of the Inner Mongolia Autonomous Region of China and the nearby border region with Russia; these areas exhibited a significant increase in high–high clustering.

4.3. The Relationships between Vegetation Dynamics and Changes in Climatic, Anthropogenic, and Hydrological Factors

As previously mentioned in Section 3.1, in order to analyze the relationship between vegetation status and climate changes, human activities and hydrological elements in the ARB, area-average climatic factors, land use area ratios, and simulated hydrological factors for 134 sub-basins—averaged annually over four phases—were used as explanatory variables to represent the environmental status of the ARB at different geographic locations and at different times. The statistical results of the explanatory variables in each sub-basin at different time periods are shown in

Table 3. Statistical results of period-average explanatory variables for four periods at the sub-basin scale.

Independent Variable		Mean	Standard Deviation	Minimum	Maximum
Climate changes	Precipitation (Pcp, mm)	565.977	100.521	289.969	800.870
	Mean temperature (Tavg, °C)	0.046	2.931	−6.253	6.051
	Maximum temperature (Tmax, °C)	6.489	2.499	0.310	12.159
	Minimum temperature (Tmin, °C)	−6.263	3.394	−12.900	0.627
Anthropogenic activities	Proportion of forest area (Forest)	0.460	0.324	0.000	0.997
	Proportion of pasture area (Pasture)	0.181	0.234	0.000	0.988
	Proportion of wetland area (Wetland)	0.049	0.081	0.000	0.542
	Proportion of crop area (Crop)	0.253	0.231	0.000	0.887
	Proportion of residential area (Residential)	0.014	0.035	0.000	0.500
	Proportion of water area (Water)	0.032	0.101	0.000	1.000
	Proportion of range area (Range)	0.011	0.065	0.000	1.000
Hydrological processes	Surface runoff (SURQ, mm)	8.941	6.423	0.000	31.384
	Groundwater flow (GWQ, mm)	2.532	3.332	0.000	32.949
	Lateral flow (LATQ, mm)	0.087	0.154	0.000	1.288
	Actual evapotranspiration (ET, mm)	33.769	7.152	14.279	76.573
	Amount of water stored in the soil profile (SW, mm)	59.386	27.181	0.000	136.027
	Snowmelt (SM, mm)	5.052	3.034	0.000	20.291

.

However, analysis of the variance inflation factor (VIF) suggested a strong multi-collinearity between these independent variables (i.e., VIF > 10). In order to obtain more stable and accurate model results, a global analysis of relationships between vegetation indexes and environmental factors in the ARB were conducted using the PLSR model [93,94]; an appropriate model suitable for cases where multicollinearity between the independent or dependent variables exists. The independent variables were filtered using variable importance in the projection (VIP) values. The PLSR model results are summarized in

Table 4. Summary of the partial least squares regression (PLSR) modeled results over the ARB.

Latent Factors (LF)	Y Variance	R2	Adjusted R2
1	0.665	0.665	0.665
2	0.097	0.762	0.761
3	0.062	0.825	0.824
4	0.008	0.833	0.832
5	0.004	0.837	0.835

, where strong (Adjusted R² >0.8) relationships between vegetation indexes and environmental factors over the ARB can be observed.

The VIP values of the independent variables are shown in

Table 5. Variable importance in the projection (VIP) values and beta coefficients of environmental factors in the PLSR model analysis.

Independent Variable	Beta Coefficient			VIP
	LAI	FVC	GPP	Model	LF 1	LF 2	LF 3	LF 4	LF 5
Pcp	0.0011	0.0002	0.4972	1.5346	1.6776	1.6435	1.5819	1.5751	1.5715
Tavg	0.0048	−0.0006	4.6085	0.7952	0.6384	0.7849	0.7568	0.7552	0.7535
Tmax	0.0078	−0.0009	6.4948	0.8564	0.7404	0.8523	0.8283	0.8274	0.8263
Tmin	0.0020	−0.0006	2.8969	0.7935	0.6432	0.7820	0.7528	0.7496	0.7481
Forest	0.6001	0.1132	237.3714	1.6577	1.7452	1.6676	1.7356	1.7357	1.7317
Pasture	−0.6216	−0.1269	−270.6097	1.4005	1.4170	1.5180	1.5098	1.5025	1.5008
Wetland	−0.6399	−0.0686	−272.2829	0.5943	0.5958	0.5673	0.5588	0.5665	0.5786
Crop	−0.2707	−0.0376	−84.2908	0.6461	0.4925	0.5442	0.5552	0.5981	0.6009
Residential	−0.9260	−0.2014	−384.5357	0.4913	0.4026	0.3803	0.4589	0.4680	0.4679
Water	−0.6076	−0.1382	−221.1102	0.5024	0.3042	0.2963	0.5232	0.5220	0.5233
Range	−0.6574	−0.1851	−237.5839	0.4536	0.4024	0.3998	0.4715	0.4755	0.4807
SURQ	−0.0052	−0.0001	−2.0949	1.2230	1.2454	1.1747	1.2252	1.2290	1.2262
GWQ	0.0020	0.0006	0.7943	0.8109	0.7918	0.7407	0.7975	0.7985	0.7967
LATQ	0.0825	0.0617	−32.9155	0.8351	0.8570	0.8210	0.8140	0.8138	0.8234
ET	0.0071	0.0008	4.4259	0.7231	0.4901	0.7509	0.7228	0.7193	0.7207
SW	0.0018	0.0003	1.1370	1.2293	1.3418	1.2865	1.2371	1.2401	1.2440
Snowmelt	−0.0006	0.0018	−0.8766	1.2590	1.3398	1.2543	1.2591	1.2552	1.2529

. Independent variables with VIP values greater than one are usually considered as relatively important in a PLSR model. For further analysis, the independent variables with VIP values greater than 0.85 (Pcp, Tmax, Forest, Pasture, Surface runoff (SURQ), Amount of water stored in the soil profile (SW), and Snowmelt (SM)) were used as inputs for the GTWR model. Regionally averaged sub-basin LAI, FVC, and GPP values were used as dependent variable inputs.

The results of the GTWR models are summarized in

Table 6. The summarized results of the ordinary least squares regression (OLS) model, temporally weighted regression (TWR) model, and geographically weighted regression (GWR) model for comparisons.

Variable	LAI				FVC				GPP
	OLS	TWR	GWR	GWTR	OLS	TWR	GWR	GWTR	OLS	TWR	GWR	GWTR
AICc	−215.4	−279.2	−745.9	−851.6	−1955.4	−2096.1	−2328.2	−2436.6	6442.5	6382.5	5665.5	5564.1
R2	0.838	0.866	0.955	0.970	0.871	0.908	0.952	0.970	0.781	0.818	0.962	0.975
Adjusted R2	0.836	0.864	0.955	0.970	0.870	0.907	0.952	0.970	0.778	0.816	0.961	0.974

NOTE: AICc: a version of Akaike information criterion (AIC) with a correction for small sample size.

. The results of the ordinary least squares regression (OLS) model, temporally weighted regression (TWR) model, and geographically weighted regression (GWR) model are also presented and were comparatively analyzed to verify the effectiveness of the GTWR model. The results indicated that the GTWR models of LAI, FVC, and GPP performed better than the other three models.

4.3.1. Regional Impacts of Climate Changes on Vegetation

As shown in Table 5, the global results of the PLSR model of the drivers prone to vegetation changes revealed the important role of Pcp in vegetation growth. Pcp is one of the factors among the independent variables that had a strong influence on LAI, FVC, and GPP (VIP > 1.5); in addition, the effects of Pcp on LAI, FVC, and GPP were all generally positive. Compared to Tmin and Tavg, Tmax and Pcp had a relatively strong effect on vegetation growth status, and the global results suggested that Tmax had a positive effect on LAI and GPP and a negative effect on FVC. Nevertheless, owning to significant spatial variability of climate changes over the study area [36,38], the relationships between vegetation and climate factors varied geographically, and thus, exhibited significant spatial heterogeneity.

The spatial distributions of the regression coefficients for Pcp with the LAI, FVC and GPP GTWR models during the four studied periods are presented in the upper part of Figure 6. Although the effects of Pcp on LAI, FVC, and GPP differed in time and space, a general tendency for the coefficients is displayed on the figures: generally lower in the northeast and relatively higher in the southwest. The areas with high coefficient values were primarily concentrated at the borders of Mongolia, Russia, and China and near the Argun River Basin and the Shilka River Basin. Areas with relatively low coefficient values were found to be mostly distributed in the northeastern part of the basin, near the border of the Sanjiang Plain in China and in Khabarovsk Krai and the Jewish Autonomous Oblast in Russia. The area of negative Pcp regression coefficient values spread over time, or the intensity of the negative effect of Pcp increased.

The spatial distributions of the regression coefficients for Tmax with the LAI, FVC, and GPP GTWR models for the four studied periods are exhibited in the lower part of Figure 6. The effects of Tmax on LAI and GPP showed similar spatial and temporal patterns and, unlike the effects of Pcp, were roughly concentrated in high-value areas in the north and low-value areas in the south of the basin. The impact of Tmax on FVC, in contrast, had a distinct area of low values in the central part of the study area.

4.3.2. Regional Impacts of Land Use Changes on the Vegetation

In addition to climatic factors, anthropogenic activities that result in land use changes, such as farming, urbanization, and reforestation, also obviously affect the vegetation dynamics in the ARB. The PLSR model results (see Table 5) revealed that at the basin scale, the proportions of forest and pasture are the most important drivers for vegetation changes. Overall, the forest area posed a positive effect on all three indexes, while an increase in pasture area most likely resulted in the decreases in all three indexes.

The upper part of Figure 7 presents the spatial distribution of the regression coefficients for forest with the LAI, FVC and GPP GTWR models for the studied four periods. From the perspective of the temporal and spatial distribution of the regression coefficients, the spatial distributions of the impacts of forest on LAI, FVC and GPP were very similar, and the spatial pattern varied slightly in the four studied periods. Regression coefficients for the forest over all sub-basins were greater than zero, and low values occurred in Huma County and Tahe County in China and Amur Oblast in Russia.

According to the spatial distributions of the regression coefficients for pasture with the LAI, FVC and GPP GTWR models, as shown in the lower part of Figure 7, the negative values of the regression coefficients were mainly concentrated near the Songnen Plain in China. In addition, areas of high pasture influence on the three vegetation indexes have decreased in size and intensity over the four studied periods. Pasture had the strongest impact on vegetation in the northeastern ARB in Period 1, i.e., 1982–1990.

4.3.3. Spatial–Temporal Heterogeneity of the Relationships between Vegetation Dynamics and Hydrological Variables

Hydrological and ecological processes are interactive: on the one hand, the water cycle strongly influences the component structure, spatial distribution, and dynamics of plant communities; on the other hand, vegetation influences hydrological processes directly through root water uptake and stomatal transpiration; while indirectly through affecting water infiltration, slope runoff and evaporation processes via the vertical canopy structure and horizontal community distribution [95]. As shown in Table 5, the PLSR models predicted that SURQ, SW and SM were more closely related to vegetation growth than the other hydrological processes considered, and had similar levels of importance. Only SW was positively correlated with all three indexes, LAI, FVC, and GPP at the global scale.

The spatial distribution of the regression coefficients for SURQ with the LAI, FVC and GPP GTWR models for the four studied periods are presented in Figure 8. Comparing the distributions of regression coefficients for SURQ with different vegetation indexes, it is clear that SURQ has the smallest range of positive regression coefficients in the FVC model. However, both the range and value of the positive coefficients in the FVC model showed an overall increasing trend during the study period.

The spatial distribution of the regression coefficients for SW with the LAI, FVC and GPP GTWR models for the four studied periods are presented in Figure 9. In the region near the Argun River Basin in the western part of the study area, the relationship between SW and LAI, FVC and GPP generally showed a roughly positive to negative trend for the four studied periods. During Period 1 (1982–1990), the region with the highest positive values of the relationship between SW and LAI/GPP appeared in the middle of the watershed, near the China–Russia border, and differed from the distribution of regression coefficients in the FVC model. The values of the regression coefficients in this region decreased with time.

Figure 10 shows the spatial distribution of the regression coefficients for SM with the LAI, FVC and GPP GTWR models for the four studied periods. In comparison with Figure 8, Figure 9 and Figure 10, the results indicate that the relationship between the hydrological factors and LAI, FVC, or GPP differed significantly in terms of the spatial pattern. The spatial distribution of negative values of the regression coefficients for SURQ with the LAI and GPP models were more similar than that with the FVC model, and the regression coefficients for SW had the similar characteristics. In contrast, the regression coefficients for SM with the LAI, FVC and GPP models were more consistent in terms of spatial distribution, with a significant increase in the range of positive values in the study area between Period 1 and Period 2 (before and after 1991) and a decrease in the area of negative values in the middle of the study area.

5. Discussion

In this study, Mann–Kendall tests and Sen’s slope were employed to examine the trends and abrupt changes in the LAI, FVC and GPP from 1982 to 2013 in the ARB. The results suggest that the LAI, FVC, and GPP exhibited strong spatial heterogeneity trends with different spatial patterns. The area demonstrating significantly increased FVC values (p < 0.05), accounting for about 62.5% of the basin area, was larger than that presenting significantly increased LAI values (7.7%) and GPP values (16.3%). The areas of LAI, FVC and GPP that exhibited significant decreasing tends (p < 0.05) accounted for 18.9, 2.0 and 10.5% of the basin area, respectively. The differences in the vegetation dynamics indicated strong spatial variations in the impacts of environmental factors on vegetation growth conditions on the one hand, and different degrees of vegetation greenness and productivity—influenced by the environment—on the other hand. However, some similarities did exist in the spatial patterns of dynamic trends of LAI, FVC and GPP from 1982 to 2013 in the ARB. For example, all three indexes showed an overall increasing trend in the Songnen Plain and a decreasing trend in the Mongolian part of the ARB. However, a much more significant increasing trend was found in the LAI and FVC over the Songnen Plain, and much more significant decreasing trends in LAI and GPP were found over the Mongolian part of the ARB. In general, changes in the environmental factors of the Songnen Plain led to an increasing trend in the greenness and productivity of the vegetation in the region.

Despite the overall increasing trend of vegetation indexes in the Songnen Plain, local Moran’s I results showed that low–low aggregation of three indexes occurred in most areas of the region during the four studied phases. For LAI and FVC, the most significant changes in local spatial autocorrelation were recognized over the Sanjiang Plain, and the low–low agglomeration in the Sanjiang Plain decreased continuously. The different characteristics and stage changes in the spatial autocorrelation in the Sanjiang Plain and Songnen Plain, both of which are the food production bases of China, may be due to different planting structures.

The PLSR results suggested that Pcp, Tmax, Forest, Pasture, SURQ, SW, and SM were the major factors affecting LAI, FVC and GPP in the entire study region. The spatial distributions of regression coefficients for the GTWR models showed that the relationships between vegetation indexes and climate factors varied with geography, and thus, exhibited significant spatial heterogeneity. In general, areas with low Pcp coefficients were found to be mostly distributed in the northeastern part of the basin, near the border of the Sanjiang Plain in China and in Khabarovsk Krai and the Jewish Autonomous Oblast in Russia. The areas with relatively high coefficient values were primarily concentrated on the borders of Mongolia, Russia, and China and near the Argun River Basin and the Shilka River Basin. A temperate humid monsoon climate predominates the eastern part of the ARB, while the Mongolian region and the Argun River Basin—as well as the Shilka River Basin in the headwaters region of the western ARB—are less influenced by the Pacific coastal monsoon and characterize a continental arid–semi-arid climate [34]. Therefore, the area near the Argun River Basin and the Shilka River Basin showed a high degree of positive influence of Pcp during all four phases, which was most likely attributable to the relatively dry climate of this region.

In addition, comparisons among the spatial distributions of the regression coefficients for Pcp with the LAI, FVC and GPP GTWR models indicated that the area and intensity of the region of high Pcp coefficient values in northeastern China decreased from Period 1 to 2 (before and after 1991). Previous studies [96,97] have reported that large areas of pasture and wetland were converted into agricultural land—especially paddy fields—in this region during this period. With significant increases in rice cultivation in this region, the extensive river runoff interception and groundwater extraction were conducted to meet the demands of irrigation [96], which probably resulted in the weakened relationship between Pcp and regional vegetation growth. In addition, the continuous conversion of pasture to paddy field in the Songnen Plain may be one of the reasons for the negative impact of pastures on LAI, FVC and GPP in the region. These results suggested that natural and anthropogenic factors jointly took effect and interacted with each other, affecting the vegetation regime of the ARB.

The variations in the spatial patterns of the regression coefficients for the hydrological variables were generally more pronounced in time than those of the climatic factors, which is most likely attributable to more dramatic spatial–temporal changes in the hydrological processes. For example, the relationship between SW and GPP shifted from positive to negative in a large area of the southwestern part of the basin during Period 2 and Period 3 (before and after 2000). In addition, a significant increase in the range of positive SM coefficient values in the study area was found between Period 1 and Period 2 (before and after 1991). To some extent, the changes represent the increasingly apparent positive effect of SM runoff on vegetation growth in the study area during the period of 1982–1999. Comparing the distributions of regression coefficients for SURQ with different vegetation indexes, the smallest range of positive regression coefficients in the FVC model indicated that the relationship between surface runoff and vegetation productivity or vegetation vertical canopy distributions may be more positive than that with vegetation horizontal cover in the ARB. The impact of a single hydrological variable on the vegetation canopy structure and canopy function may vary in both direction (positive or negative) and intensity.

6. Conclusions

This study analyzed the spatial–temporal variations and spatial autocorrelation of vegetation growth conditions in the ARB at the pixel scale and investigated the spatial–temporal heterogeneity of the relationships between vegetation indexes and climatic, anthropogenic and hydrological factors at the sub-basin scale, using GTWR models from 1982 to 2013.

The trends of LAI, FVC and GPP in the ARB from 1982 to 2013 were quite different, but only the variation of area-average annual FVC proved to be significant with an increase rate of 0.0004/year (p < 0.05). All the annual vegetation indexes underwent a rapid change around the year 2001. Additionally, all three indexes showed an overall increasing trend in the Songnen Plain and a decreasing trend in the Mongolian part of the ARB during 1982–2013. The local Moran’s I results of LAI, FVC, and GPP suggested that similar decreasing trends in low–low agglomeration area of LAI and FVC on the Sanjiang Plain occurred during the study period.

According to the results of the GTWR models, the relationships between vegetation indexes and driving factors exhibited significant spatial–temporal heterogeneity. For example, a general tendency of the coefficients for Pcp existed: generally lower in the northeast and relatively higher in the southwest. The area near the Argun River Basin and the Shilka River Basin showed a high degree of positive influence on Pcp during all four phases, which was most likely attributed to the relatively dry climate in this region. Additionally, the negative values of regression coefficients for pasture were mainly concentrated near the Songnen Plain. The continuous conversion of pasture to paddy field may be one of the reasons for the negative impact of pastures on the vegetation dynamics in this region. SURQ had the smallest range of positive regression coefficients in the model for FVC among the three vegetation indexes, indicating that the relationship between surface runoff and vegetation productivity or vegetation vertical canopy distributions may be more positive than that with vegetation horizontal cover. Additionally, a significant increase in the range of positive regression coefficients for SM between Period 1 and Period 2 may suggest the increasingly apparent positive effect of SM runoff on vegetation growth in the study area during the period of 1982–1999.

Natural and anthropogenic factors jointly took effect and interacted with each other to affect the vegetated regime of the region. The decrease in the impact of precipitation on vegetation growth in the Songnen Plain was determined as having commenced around 1991, which was most likely attributable to dramatic changes in water use styles induced by local land use changes, and this corresponded to the negative correlation between pasture areas and vegetation indexes during the same period.