Spatial and Temporal Characteristics of Vegetation NDVI Changes and the Driving Forces in Mongolia during 1982–2015

: As a result of the unique geographical characteristics, pastoral lifestyle, and economic conditions in Mongolia, its fragile natural ecosystems are highly sensitive to climate change and human activities. The normalized difference vegetation index (NDVI) was employed in this study as an indicator of the growth status of vegetation. The Sen’s slope, Mann–Kendall test, and geographical detector modelling methods were used to assess the spatial and temporal changes of the NDVI in response to variations in natural conditions and human activities in Mongolia from 1982 to 2015. The corresponding individual and interactive driving forces, and the optimal range for the maximum NDVI value of vegetation distribution were also quantified. The area in which vegetation was degraded was roughly equal to the area of increase, but different vegetation types behaved differently. The desert steppe and the Gobi Desert both in arid regions have degraded significantly, whereas the meadow steppe and alpine steppe showed a significant upward trend. Precipitation can satisfactorily account for vegetation distribution. Changes of livestock quantity was the dominant factor influencing the changes of most vegetation types. The interactions of topographic factors and climate factors have significant effects on vegetation growth. In the region of annual precipitation between 331 mm and 596 mm, forest vegetation type and pine sandy soil type were found to be most suitable for the growth of vegetation in Mongolia. The findings of this study can help us to understand the appropriate range or type of environmental factors affecting vegetation growth in Mongolia, based on which we can apply appropriate interventions to effectively mitigate the impact of environmental changes on vegetation.


Introduction
As an important component of the surface ecosystem, land vegetation can fundamentally regulate the energy balance, and water and biogeochemical cycles of the Earth's surface through photosynthesis, respiration, transpiration, and surface albedo [1,2]. Furthermore, land vegetation is largely dependent on and sensitive to a variety of factors of natural environment and human activities, reflecting the effects of climate change and human activities over a short period of time [3,4]. Vegetation also affects the economic structure and economic development of a country or region, especially in arid and semi-arid areas, where agricultural and livestock production are the main economic activities [5]. Therefore, investigating the relationship between land cover change with respect to climate change and human activities is important for assessing economic development and understanding ecosystems.
With the continuous development of earth observation system technology and the continuous enrichment of statistical models, the research content and methods of the relationships between vegetation, climate change, and human activities have become increasingly diverse. The normalized density vegetation index (NDVI), derived from infrared channel and near-infrared channel remote sensing data, is a good indicator of vegetation growth status and spatial distribution density of vegetation, and is linearly related to vegetation distribution density [6,7]. A large number of ecological studies have been carried out related to vegetation change and its influencing factors using NDVI and statistical models. Temperature and precipitation, as important factors in the natural environment, are often used as research priorities to assess their impact on vegetation. However, the dominant factors contributing to NDVI in different regions vary considerably. The variety in vegetation and its connection to climate change have been investigated by Du et al. [5] and Zheng et al. [8] in China. These studies identified precipitation as the key climatic factor governing variation in NDVI. However, Guo et al. [9] found that temperature was the dominant factor on vegetation (NDVI) growth during the growing season in permafrost zone of Northeast China. Baniya et al. [10] also found that vegetation (NDVI) change in Nepal was more sensitive to temperature than precipitation change. Conversely, Peng et al. [11] found that the soil type was the governing natural factor on NDVI changes in Sichuan, western China. These studies indicate that the diversity of geographical environment determines the dominant factors affecting vegetation growth and a limited number of factors can influence vegetation growth.
For studies such as those mentioned above, correlation analysis, partial least squares regression, and multivariate regression analysis or their disguised methods (e.g., the RESTREND method [12]) are the most frequently used analysis techniques. However, because of the complex process of the land vegetation growth response to driving factors, the inflexible statistical linear models may not accurately describe the internal relationships between the two variables [12]. Spatial heterogeneity is a pervasive feature of ecological and geographical phenomena, which has been referred to as the Second Law of Geography [13][14][15]. Techniques for spatial analysis need to consider spatial heterogeneity. Traditional research methods tend to unify the relationships between independent variables and dependent variables, which is not in accordance with the Second Law of Geography [16]. Based on the theory of spatial (global) stratified heterogeneity, Wang [16] proposed the geographical detector model (GDM), which was designed to measure the spatially stratified heterogeneity of a response variable and reveal the impact of driving factors. The GDM, which is independent of any linear hypothesis, can detect the nonlinear relationships between two variables from the perspective of spatial heterogeneity and stratification.
Due to the inland environment in northern Central Asia, Mongolia is located in the tail-end of monsoons of the Pacific Ocean and the Arctic Ocean. Most parts of this county are arid and semi-arid regions with fragile ecological environments and are sensitive to climate change and human activities. Over the past 20 years, the temperature increase rate has been about three times the rate of global average temperature rise (0.06 °C/a) and precipitation has shown a decreasing trend with a rate of change of −1.4 mm/a in Mongolia [17,18]. Historically, herding has long been the dominant agricultural activity in Mongolia. In the early 1990s, Mongolia experienced a social institutional change, and since then, livestock quantities have grown at an unprecedented rate, especially in the north-central part of Mongolia [19,20]. According to the National Statistical Office of Mongolia, the number of livestock has increased nearly 2.5 times (from 25,000 million to 66,000 million) over the last 30 years. Therefore, the spatial and temporal distribution of climate and economic activities in Mongolia is undergoing major changes leading to high spatial heterogeneity. In addition, a few studies have shown that Mongolian vegetation is particularly vulnerable to climate change and human activity [21][22][23].
To date, numerous studies have investigated the driving factors of the vegetation change in Mongolia at different spatial and temporal scales; however, some crucial issues have still not been fully addressed. Firstly, most analytical methods on the factors influencing NDVI change were limited to simple linear regression or correlation analyses [20,21,24,25]. Given the complex relationships between geographic variables described above, these linear hypotheses usually lead to biased results. Second, previous research commonly explored the relationships between vegetation growth and climatic factors from a time series, ignoring their spatial differences [26][27][28][29][30]. Thirdly, previous studies have not quantitatively assessed the interaction between two or multiple factors, which are commonly used to quantitatively check whether two or multiple environmental determinants work independently or not. It is extremely important therefore to conduct a comprehensive evaluation of the effects of environmental factors on vegetation growth.
In this study, we first conducted a spatiotemporal analysis of vegetation changes in Mongolia from 1982 to 2015. Then, we identified the main environmental factors (natural factors and human activity factor) and their relative roles in vegetation distribution. In addition, we also measured the driving force of environmental factors on vegetation changes. We determined the interaction between factors and the optimal range of each factor beneficial to vegetation growth in Mongolia. Study of vegetation growth and change of vegetation and its relationship with environmental factors will help pasture managers to select appropriate grazing activities from the perspective of sustainable pasture use. This research can also provide a reference for conservation managers. We would expect that targeting of the dominant factors could reduce the investment cost of vegetation restoration projects.

Study Area
Landlocked Mongolia is situated in the north-central Asian mainland, deep within the inland far from any ocean (41.5°-52.1°N, 87.70°-119.9°E), and is basically a plateau with an average elevation of about 1580 m above sea level (Figure 1a). The special geographical environment of Mongolia leads to clear continental climate characteristics: a cold long winter and hot summer, with large annual and diurnal temperature ranges. Mongolia has three belts of varying basic vegetation type: forest, steppe, and Gobi Desert from north to south (i.e., from mountains to the basins and plains). The spatial distribution of vegetation coverage is roughly parallel to the latitude (Figure 1b). The steppe is the main land cover type of Mongolia, accounting for 66% of the total area of the country. The spatial distribution of forests is concentrated, mostly in the mountains of northern Mongolia, and the steppes are mainly found in central Mongolia and distributed along latitudes. The Gobi Desert area is mainly distributed in the southwest of Mongolia. Vegetation is sparse and rare in this area, consisting mainly of herb wormwood and halophytic vegetation. For convenience, we collectively refer to all vegetation in this area as Gobi Desert. With a total population of 3.2 million (January 2019, National Statistics Office of Mongolia: http://www.1212.mn/) and land area of 1.56 × 10 4 km 2 (Figure 1b), Mongolia has one of the lowest average population densities in the world. Economic activity in Mongolia has long been based on herding. According to the statistical data from the National Statistics Office of Mongolia in 2019, about 25% of Mongolia's population drive their incomes from breeding livestock. Most herders in Mongolia follow a pattern of nomadic or semi-nomadic pastoralism.

Data Sources
In this study, we used NDVI, climatic, topographical, livestock quantity, soil type, and vegetation type data (Table 1).  [24,32,33]. The spatial and temporal resolution of this dataset is 1/12° and bimonthly frequency, respectively. Different NDVI extraction methods can obtain different vegetation conditions. The widely used NDVI composite methods include average value composite (AVC) and maximum value composite (MVC) methods. The average value composite method integrates all vegetation states within a year or during the growing season and can more comprehensively reflect the vegetation year information. However, the average value composite method is affected by many factors, such as growth period, cloud, and atmospheric conditions. In this study, we used the maximum value composite (MVC) method to preprocess the NDVI3g dataset from bimonthly to annual. The MVC method was used to extract the maximum NDVI value in the multi-phase data at each spatial location as the real NDVI value [34]. MVC also has some limitations affected by the type of vegetation and growth period. For example, with the sudden increase of precipitation in summer, the rainfed grasslands of southern Mongolia may have extremely high NDVI in a short time, and then its NDVI value is almost close to zero for the rest of the year. Taking the highest NDVI that occurred in this year as the final NDVI value ignores the seasonal changes of vegetation. Considering that this research is devoted to investigating the trend and driving force of the interannual change of NDVI, the NDVI obtained by the MVC method is more stable in the performance of the interannual change of vegetation. The spatial distribution of the annual average NDVI of Mongolia from 1982 to 2015 is shown in Figure 1c.

Climatic Data
The climatic data used in this study include precipitation, air temperature, wind speed, snow depth, and specific humidity (Figure 2a-e). All the climate data were procured from Global Land Data Assimilation System (GLDAS) Catchment Land Surface Model L4 V2.0 at daily, 0.25° × 0.25° resolution, which belong to the Goddard Earth Sciences Data and Information Services Center (GES DISC) [35]. The CLM (Catchment Land Surface Model) data set has been used in a variety of application studies such as soil moisture [36], permafrost [37], and groundwater storage changes [38].
In addition, we collected extra data from 70 meteorological stations in Mongolia (Figure 1d) to evaluate the data quality of the GLDAS CLM. We did not use the meteorological data directly in our study because of the insufficient number of available meteorological stations in Mongolia and also because the spatial interpolation results cannot accurately represent the real local meteorological state. The GLDAS CLM data are assimilated with a combination of model and observation-based forcing datasets, which is more accurate than direct interpolation. To further verify the reliability of the GLDAS CLM, we performed a simple agreement assessment on the precipitation (Total precipitation rate) and temperature (Air temperature) data of CLM and meteorological stations. The assessment methods were coefficient of determination (R 2 ), the Pearson correlation coefficient (r), and root mean square error (RMSE). To evaluate the overall agreement between CLM data and meteorological station data, we compared the two datasets for all pixels, which included all the meteorological stations (70 stations) in space and time (34 years from 1982 to 2015). The assessment results indicate that the R 2 values of precipitation and temperature were 0.7 and 0.8, respectively, the r values were 0.83 and 0.89, respectively, and the RMSE values were 68 mm and 1.6 °C. The results indicate that there is a high agreement between CLM data and meteorological station data. These data were used as input data for our model. In Mongolia, annual average temperature ranges from 12.5 °C in the meadow steppe to 10.29 °C in the Gobi Desert and annual average precipitation varies from 50.8 mm in Gobi Desert to 596.9 mm in forest [17,18].

Topography
The topographic data used in this study included elevation, slope degree, slope aspect, and curvature (Figure 2f-i). All topographic data were derived from the digital elevation model (DEM). The DEM data we selected in this study were from Shuttle Radar Topography Mission (SRTM) Version 4 with 90 m spatial resolution. This is now the highest quality SRTM dataset available with a horizontal standard error of 1 m and vertical standard error of 0 m [39].

Livestock
Herding of livestock is the primary economic activity of the Mongols. About three-quarters of the land in Mongolia is made up of pastures that provide support for an enormous quantity of grazing livestock. Due to the small population density of Mongolia, the impacts of industry, agriculture, and urban construction on vegetation are limited [24]. Grazing activities are the primary and most direct factors affecting vegetation change. Therefore, in this study, we hypothesized that the quantity of livestock could be used as a human activity factor to assess the influence on vegetation. Livestock quantity data were obtained from the National Statistics Office of Mongolia ( Figure 2j). The livestock quantity was collected at city level. In this study, livestock refers to horses, cattle, camels, sheep or goats, which are also the five largest livestock species grazing in Mongolia.

Sen's Slope
Sen et al. [41] proposed an estimator to analyze the slope trend of a time series data. Sen's slope is a non-parametric method requiring the test data independent and insensitive to missing and abnormal data. The Sen's slope is computed as follows: where xi and xj are the values at times i and j, respectively, and 1 ≤ i < j ≤ n, n refers to the time series data length (i.e., 34 years from 1982 to 2015) in the current study. The slope sign reflects data trend, indicating the rate of variation of the time series data. Slope > 0 means the upward trend and the Slope < 0 indicates a downward trend. Since this method lacks statistical significance testing for the trend, we assessed it using the Mann-Kendall test method, which also a non-parametric statistical test.

Mann-Kendall Test
Mann et al. [42] and Kendall et al. [43] introduced the Mann-Kendall (M-K) test for evaluating the significance of trends of time series data. As a non-parametric statistical test, it is also not affected by outliers. The M-K test is often used to quantify the significance of hydrometeorological time series trends [11,[44][45][46][47]. The M-K test is given below: where n refers to the time series data length, and xi and xj are the data values in the time series data at times i and j, respectively. When n > 10, the statistic S approximately equals the standard normal test statistic (Z) value, which can be used to test the trend as following: where n refers to the time series data length, m is the number of repeated datasets in the time series data and t denotes the repeated data values in the ith group. When |Z| < Zα, the null hypothesis can be rejected, meaning the existence of a significant trend in this data. When |Z| > Zα, the null hypothesis of trend will be accepted. In this study, α = 0.01 and α = 0.05 are defined as the given significance levels, |Z| − α/2 is equal to 2.576 and 1.96, respectively.

Geographical Detector Model
The geographical detector model (GDM) proposed by Wang et al. [16] is a novel tool that not only can quantify the relationship between two variables Y (response variable or geographical phenomenon) and X (explanatory variable or environmental factor), but also can investigate the interaction relationship between two explanatory variables (X1 and X2) to the explanatory variable Y, according to spatial distributions of these data, without assumption of linearity of the association [16,48]. The core idea of GDM is that if an explanatory variable (Xi, i = 1, 2, ...) contributes to a response variable (Y), the spatial distribution between the two is of similarity, and this similarity or spatial association can be measured by the power of determinant (q) value [16] or q-statistic [48]: (6) where N and σ 2 refer to the number and the variance of Y values; Y is composed of L strata (h = 1, 2, ..., L); Nh and σ 2 denote the number and the variance of Y values in the hth stratum. The value of q is strictly within [0,1]. The q value equals 0, there is no association between X and Y, whereas if q equals 1, Y is completely determined by X. The geographical detector is composed of four subdetectors: factor detector, interaction detector, risk detector, and ecological detector. For more information about geographical detector model, readers can refer to the website of Geodetector: http://geodetector.cn/.
The factor detector is used to quantify the driving force of each environmental factor X on response variable Y. If X is numerical, the first step is to discretize X into several strata depending on their values. There are many discretization methods, such as natural breaks method and equal interval method. The natural breaks method produces discrete intervals that have the characteristics of the smallest in-class variance and biggest variance between classes, which designed by Jenks et al. [49]. The equal interval method equally divides data into specified quantities. The discretization method used in this study is a classification method based on the principle of maximum q value, which can uncover more correlative information between environmental factors and geographical phenomenon [16].
The interaction detector reveals whether the environmental factors X1 and X2 (or more) have an interactive influence on a response variable Y. This relationship evaluation method involves comparing q(X1), q(X2), and q(X1 ∩ X2) including weakening, independent, and enhanced relationships (Table 2). Meanwhile, interaction detectors can measure whether two factors weaken or enhance each other, or whether they are independent of vegetation growth. Therefore, the interaction results not only indicate the contribution of the two factors to the growth of vegetation, but also reveal whether there is co-linearity (correlation) between the two factors. For example, when the interaction result is q(X1 ∩ X2) < Min(q(X1), q(X2)), indicating that the two factors shoe a nonlinear weakened effect and a strong correlation between the them. Similarly, when the interaction result is q(X1 ∩ X2) > q(X1) + q(X2) or q(X1 ∩ X2) = q(X1) + q(X2), indicating that the two factors show a nonlinear enhanced effect or independent effect and there is no correlation between the two factors.

Description
Interaction The risk detector is used to identify the favorable and unfavorable strata of environmental factor to the response variable. It is implemented by measuring the significance of difference between two strata of the environmental factors using a t-test. The t-test is calculated as follows: where Y, σ 2 , and N refer to the average, variance, and number of stratum of data Y. The null hypothesis means that there is no difference between strata 1 and 2. If α, we accept the null hypothesis and we need to examine the underlying natural mechanisms.
The ecological detector is used to judge whether one environmental factor is more significant than another. For example, in this study, ecological detector can compare the effect of precipitation and temperature on the vegetation growth, and we can identify which effect is more significant. The test method can be found where F denotes the F-test result value, and Nx1 and Nx2 refer to the number of environmental factors X1 and X2, respectively. The σ 2 denotes the variance of the impact factor. The significance level used herein is 0.05.

Dynamic Variation of NDVI
The proportion of NDVI change with an increasing trend (19%) from 1982 to 2015 approximately equaled the proportion with a decreasing trend (18%), although regional differences were found ( Figure 3). The areas with significant variation in vegetation change were mainly located in eastern regions (Dornod and Sukhbaatar province), southern regions (Govi-Altai, south of Bayankhongor and west of Umnugovi province), western regions (Bayan-Ulgii province), and central regions (Arkhangai province) of Mongolia. Among these regions, the eastern and western regions showed a clear greening tendency. In contrast, the southern and central regions showed the greatest decrease.  Figure 4. In this part of the analysis, we only counted those pixels that were statistically significant. For forest, pixels exhibiting an increasing trend and those with a decreasing trend occupied 17% and 15.4% of the total forest area, respectively, implying that these two trends each occurred over a similar areal extent. The highly significant (p < 0.01) and significant (0.01 < p < 0.05) trend of increasing and decreasing were also similar, indicating that the forest vegetation type in Mongolia showed no unified changing trend between 1982 and 2015.
For meadow steppe, areas with an increasing trend and those with a decreasing trend occupied 31.4% and 3% of the total steppe area, respectively, indicating that the area of Mongolian meadow steppe under an increasing trend was much higher than that under a decreasing trend. The area with a highly significant increasing trend was higher than that with a significant increasing trend. However, the area with a highly significant decreasing trend was lower than that with a significant decreasing trend. These findings suggest that meadow steppe in Mongolia went through a general positive development during 1982 to 2015.
For typical steppe, the proportion of pixels exhibiting an increasing trend and those with a decreasing trend were 17.1% and 9.9% of the whole typical steppe area, respectively, suggesting that the area under an increasing trend was higher than that under a decreasing trend during 1982-2015. The area with a highly significant increasing trend was lower than that with a significant increasing trend. In contrast, the area with a highly significant decreasing trend was higher than that with a significant decreasing trend. These findings imply that typical steppe in Mongolia has shown increased vegetation greening during this period. For desert steppe in Mongolia, the proportion of pixels exhibiting an increasing trend and decreasing trend were 5.6% and 12.4%, respectively, suggesting that the area under an increasing trend was much lower than that under a decreasing trend. Both the areas of highly significant increasing and decreasing trends were lower than those with a significant level, which implies that desert steppe in Mongolia showed vegetation degradation during this period.
For alpine steppe, the proportion of pixels exhibiting an increasing trend and decreasing trend were 26.3% and 5.8%, respectively, demonstrating that the area under an increasing trend was much higher than that under a decreasing trend. The areas of highly significant increasing and decreasing trends were lower than those with a significant trend, which implies that alpine steppe in Mongolia experienced increased greening during this period.
For Gobi Desert in Mongolia, the proportion of pixels exhibiting an increasing trend and decreasing trend were 8% and 35.1%, respectively, indicating that the area under an increasing trend was considerably lower than that under a decreasing trend during 1982-2015. The areas with highly significant increasing and decreasing trends were higher than those with a significant level, suggesting that Gobi Desert in Mongolia experienced a severe degradation process between 1982 and 2015. In order to further understand the NDVI changes in Mongolia from 1982 to 2015, we focused on the relationships between the NDVI values at the beginning of the change and the variation within the 34-year period (Variation = Sen's slope × 34). Figure 5a shows the scatter diagram of significant points (p < 0.05) calculated by the Sen's slope method. The X-axis of Figure 5a reveals that the initial value range is between 0.05 and 0.95. However, based on the statistical results of the histogram of the initial value, the initial value of the changed region is mainly distributed within two levels: 0.05-0.2 and 0.3-0.5. From a spatial perspective, this indicates that the area of southern Mongolia containing Gobi Desert areas underwent major vegetation changes from 1982 to 2015. From the Y-axis, the variation is mainly distributed between −0.2 and 0.2, and the variation less than 0 (i.e., vegetation degradation) is mainly concentrated at −0.1 to 0. Combined with the specific spatial location, although the southern Gobi Desert region mainly reduced over the 34-year period, the spatial reduction was small and the amount of vegetation degradation was also small. However, from an ecosystem perspective, this degradation has a relatively large impact on the southern grassland ecosystem. In Figure 5a, the variation count histogram reveals that, overall, the amount of vegetation increase was less than the amount of reduction. This contrasts with the areal analysis shown in Figure  3, in which the area of vegetation increase (15% of pixels) was larger than the area of reduction (14% of pixels). The distribution of reduction was wider (interval of 0-0.2), whereas the increase was more concentrated. The NDVI variation of different vegetation types is shown in Figure 5b. Meadow steppe, typical steppe, and alpine steppe vegetation types showed an increase, whereas the amount of variation in meadow steppe and desert steppe showed a decrease.

Trends of the Environmental Factors
The Sen's slope and M-K test were applied to the 34-year period for the environmental factors including five climate factors and one human activity factor: precipitation, air temperature, wind speed, snow depth, specific humidity, and livestock quantity. The output of spatial distribution of these environmental factors' trends and the corresponding areas were analyzed ( Figure 6). As demonstrated in Figure 6a, precipitation showed a high spatial heterogeneity of different trends across Mongolia. The slight increasing precipitation appeared in the north (Khuvsgul province) and south (Uvurkhangai province) of Mongolia with a rate of change of approximately 2 mm/a. The distinct decrease zones of precipitation were detected at the middle area (Ulaanbaatar, Khentii and Tuv province) of Mongolia, with rate of change of −4 mm/a. In general, most areas of Mongolia showed a decreasing trend in precipitation during 1982-2015 (Figure 6a).
In contrast, the temperature in Mongolia presented the opposite areal trend. As shown in Figure  6b, the temperature across Mongolia followed an observably increasing trend during 1982-2015, with no areas of decline. The highest temperature rise rate reached 0.06 °C/a, which is consistent with the results of previous studies [17]. The highest temperature increases were found in the southeast and northwest of Mongolia and the lowest temperature increases were found in the north and northeast boundary of Russia and Mongolia. Overall, Mongolia experienced a rapid and accelerated warming trend from 1982 to 2015.
Trends in wind speed factor also exhibited high spatial variability in Mongolia. As shown in Figure 6c, the largest increase and decrease were found at the northwest (Khovd, southern Zavkhan, and northern Govi-Altai province) and center (Uvurkhangai, Umnugovi, and Dornogovi province) of Mongolia, respectively, with an approximate trend of 0.015 m/s per year.
The snow depth of Mongolia mainly showed a positive trend (Figure 6d), with an increase rate of 7.6 mm/a. Few areas with a negative trend occurred across southeast and west Mongolia.
There was clear spatial trend distribution of specific humidity across Mongolia (Figure 6e). The specific humidity in the east of Mongolia generally decreased, whereas the specific humidity in the west area generally increased. The extreme values of these two trends were −20 and 13 mg/kg per year, respectively.
From 1982 to 2015, Mongolia experienced a considerable growth in livestock quantity. According to the National Statistical Office of Mongolia, livestock quantity increased substantially from 24.76 to 55.97 million during 1982-2015 and reached up to 66.21 million by 2018. From the perspective of spatial trend (Figure 6f), an overwhelming majority of city level showed a rapid increase in livestock quantity, and those with a relatively high positive trend were mainly found in the north-central and northeastern parts of Mongolia, with the increase rate of 15 thousand/a.

Single Factor Influence Detection
The factor detector of GDM was used to reveal the driving force of each environmental factor on response variable NDVI. The higher the q value obtained from factor detector, the stronger the contribution of the factor to response variable. Moreover, the factor with maximum q value was defined as dominant. In this analysis, the response variable was taken as the annual average value of NDVI during 1982-2015. The climate factor and livestock quantity factor herein were also taken as annual averages for 1982-2015. All of the environmental factors and response variables are shown in Figure 6a. The results of the calculation of the q value of each environmental factor to NDVI are shown in Table 3. All factors were ranked in descending order by their dominant power as follows: Prec > Vegett > Soilt > Snowd > Temp > Winds > Slopd > Livstq > Specfh > Elev > Curv > Slopa (first row in Table 3).
Among the different environmental factors, the q values of precipitation, vegetation type, and soil type were the three largest factors (Table 3). Therefore, precipitation is the dominant factor affecting vegetation distribution in Mongolia. The q values of snow depth, temperature, and wind speed were 0.49, 0.46, and 0.45, respectively, and these together accounted for more than 40% of the NDVI distribution. However, the slope degree, specific humidity, elevation, curvature, and slope aspect affect the spatial distribution of hydrothermal conditions and thus affect the growth of vegetation. The influence of each factor is indirect and all q values for the NDVI distribution were below 20%, demonstrating that the factors of slope degree, specific humidity, elevation, curvature, and slope aspect had only a minimal influence on the NDVI distribution. Notably, the q value of human activity factor (livestock quantity) was 0.21, which indicates that grazing has little effect on the distribution of vegetation when considering Mongolia in its entirety. All driving factors of vegetation distribution test passed the significance test (p < 0.05)

Combined Influences Detection
The interaction between pairs of environmental factors (symbolized by ∩) was analyzed and the results are shown in Table 4, including comparison of their interactive q value and individual q value to NDVI. The top five interactive q values decreased in the following order: Prec ∩ Vege > Prec ∩ Soilt > Vegett ∩ Soilt > Prec ∩ Winds > Prec ∩ Slopd. This indicates that the interactions between meteorological factors (precipitation or wind speed), vegetation type, soil type, and topography had the greatest impact on the NDVI value. Two factors exhibiting a nonlinear enhancement after the interaction means that the q value after interaction is greater than the sum of the individual q values before interaction (q(X1 ∩ X2) > q(X1) + q(X2)). In this study, factors that showed nonlinear enhancement after interaction included Snowd ∩ Elev, Spec ∩ Slopd, Winds ∩ Elev, Spec ∩ Elev, Elev ∩ Slopd, Winds ∩ Curv, Spec ∩ Curv and Spec ∩ Slopa. This indicates that the influence of topographic factors (Elev, Slopd, Slopa, and Curv) as indirect factors on vegetation growth are mainly reflected in the interaction with other factors. For example, when interacting with meteorological factors, topographic factors can markedly enhance their single effects on the NDVI value. Naturally, all interactive q values of pairs of environmental factors were larger than any q values of an individual factor, indicating that the effects of environmental factors on the distribution of NDVI are not independent but interactive. Compared with these factors that enhanced each other, the correlation between these nonlinear enhancement factors is relatively low (Table 4). Note: "↑" and "↑↑" refer to X1 ∩ X2 enhancing each other and X1 ∩ X2 nonlinear enhancement, respectively.

Optimal Range of Factors for NDVI Detection
Using the risk detector of GDM, the optimal range results for the main environmental factors are shown in Table 5. We assumed the limits of environmental factors with the biggest average value of NDVI is the optimal range for vegetation growth. There was a range of different relationships between different factors and the vegetation distribution. For example, there was a positive correlation between precipitation and NDVI value, indicating that the annual average NDVI value increased with the increasing precipitation. In this study, when the annual precipitation was between 331 mm and 596 mm, the annual average NDVI value reached maximum. In contrast, the NDVI increased with decreasing temperature and the maximum NDVI value appeared in the first stratum (stratum for minimum temperature) of temperature.
Only monotonic relationships were found between precipitation, temperature, and vegetation growth; other factors were non-monotonic. For instance, the second stratum of wind speed (2.74-3.27 m/s) contained the maximum NDVI, implying that this stratum of wind speed promoted widespread vegetation growth in Mongolia. As for snow depth, the maximum annual average NDVI value appeared in the stratum between 92 mm and 94 mm, implying that this stratum was the optimum for vegetation growth. The maximum annual NDVI value of 0.49 roughly corresponded to an altitude between 535 m and 974 m. We also found that NDVI values can reach the peak value when the stratum of specific humidity is between 4.51 mg/kg and 5.25 mg/kg. As for slope degree, slope aspect, and curvature, the optimal range of these factors were 4.1-23.8°, 330-360°, and −0.62-0.08, for which the NDVI can reach 0.53, 0.43, and 0.49, respectively. For the vegetation types, the highest NDVI value was found for the forest, followed by meadow steppe. For the soil types, the highest NDVI value was found for the pine sandy soil, followed by mountain peat soil and grey forest soil, in which all the average NDVI values exceeded 0.70. For livestock quantity, the NDVI peak occurred in the stratum between 128,000 and 132,000. These findings could play a guiding role in restoring vegetation engineering in arid land. The detection results of the dominant range of environmental factors in this study reveal the dominant habitat range for vegetation growth, which can help land managers in planning a suitable range of vegetation growth environment. It also plays the role of planting grass or cash crops as much as possible with the minimum human intervention. In short, the detection results of the dominant range of environmental factors can provide a basic environmental framework for vegetation growth.

. Significant Differences between Factors
The significant differences between the effects of environmental factors on vegetation growth from the ecological detector in GDM are shown in Table 6. The topographical elevation, slope degree, and slope aspect all had the same influence on vegetation growth. In contrast, for the slope degree and snow depth factor, only wind speed and temperature factors had no significant effect on the spatial distribution of vegetation, respectively. There were significant differences between all of the other environmental factors. Table 6. Statistical significance of the detection factors (95% confidence level).

Factor Prec Soilt Snowd Temp Winds Vegett Slopd Specfh Elev Cuve Slopa
Note: "Y" and "N" refer to whether there is Yes or No significant difference between two factors.

Single Factor Influence Detection
Using the annual average NDVI value as the response variable, the annual average climate factors and the topographical factors, as well as the livestock quantity as environmental variables to explore their internal relationships is inclined to reveal the physical mechanism between them. Therefore, we used the Sen's slope method to obtain the time series change trend of NDVI and environmental factors, and then used the GDM to detect the relationships between their changes. All these slope variables are shown in Figures 3 and 6. Note that only six factors changed over time, including five climatic factors and the livestock quantity factor. In order to ensure a complete contrast, the invariant factors (i.e., topographical factors, vegetation type, and soil type) were also included in the model calculations.
The results of the factor detector changes in GDM are shown in Figure 7a (blue curve) and Table  3 (the second row). All factors were ranked in descending order by their dominant power: Elev (0.2026) > Temp (0.1232) > Livstq (0.122) > Specfh (0.1091) > Vegett (0.0984) > Soilt (0.098) > Prec (0.0473) > Winds (0.0318) > Slopa (0.0206) > Slopd (0.0135) > Curv (0.0131) > Snowd (0.0074), indicating that vegetation changes during 1982-2015 were mainly due to changes of elevation and temperature, rather than the precipitation and vegetation type which are the dominant controlling factors on vegetation distribution. Furthermore, changes in livestock quantity ranked third in terms of influence on vegetation changes, but eighth for its influence on vegetation distribution. For Mongolia as a whole, the effects of these factors on the vegetation change were minimal (all the q values were less than 0.21), indicating that the changes of environmental factors do not match the changes of vegetation very well in terms of spatial distribution across Mongolia. This is because the physical mechanism is complicated and spatial heterogeneity is ignored when regarding the entirety of Mongolia as a single research area. This may lead to misunderstanding of the analysis results.
However, our study uses vegetation types as the basic unit to explore the spatial and temporal characteristics of NDVI changes and the corresponding drivers, which are more in line with the natural geographical law. Therefore, our statistical results can approximately reveal the relative influence power between factors on vegetation changes and could be used as a reference with which to compare the contributions of different factors.

Vegetation Change Drivers under Different Vegetation Types
Since the analysis of the whole of Mongolia as a research area ignores its spatial heterogeneity, we used a partitioning approach to separately analyze which environmental factors drive vegetation changes and thus provide a more accurate assessment of the rate of change between vegetation and environmental factors. In this section, we detected factors based on the six different vegetation types (Figure 2k). The results of factor detection are shown in Figure 7b.
For forest, the maximum q value was for livestock quantity (q = 0.5), followed by elevation (q = 0.22) and wind speed (q = 0.17). Similarly, for meadow steppe, the dominant factor was livestock quantity (q = 0.53), followed by elevation (q = 0.37) and precipitation (q = 0.31). For the typical steppe, the dominant factor was elevation (q = 0.26) and the secondary driver was soil type (q = 0.17). For the desert steppe, the main influencing factor was also the livestock quantity (q = 0.27). The contribution power of elevation and specific humidity were 0.19 and 0.13, respectively. For Alpine steppe, the dominant factor was slope aspect, for which the q value could reach 0.46. In addition, other factors also had high driving force values, such as elevation (q = 0.43), temperature (q = 0.39), and slope degree (q = 0.33), indicating that these factors have major impact on vegetation changes. For Gobi Desert, as shown in Figure 7b, the dominant factor was again livestock quantity (q = 0.44). However, the remaining factors such as precipitation (q = 0.15) and wind speed (q = 0.14) showed a weak influence on vegetation changes. Overall, based on the analysis of different vegetation types, the dominant factors that affect vegetation change were livestock quantity and topographical factor.

Discussion
The analysis of the driving factors and trends of vegetation growth and vegetation changes derived from satellite observations are of great importance for evaluating the value of grassland ecosystems in Mongolia. However, a major challenge remaining for research into the driving force of vegetation growth and changes in Mongolia is that analytical models or methods fail to take the spatial heterogeneity of variables into account. In this study, we used Sen's slope and the M-K test to allow trend analysis of vegetation changes. These two pixel-based analysis methods can be used to delve into spatial heterogeneity information of vegetation changes. Then, the GDM based on spatial stratified heterogeneity [48], was utilized to investigate the driving force of vegetation growth and changes.

The Trend of Vegetation Changes
Bao et al. [21] used the GIMMS NDVI to study the vegetation dynamics in Mongolia, and found that different vegetation types had different variation trends, but only the forest and alpine steppe exhibited a significant increasing trend and the meadow steppe, typical steppe, and desert steppe exhibited no significant trend. The trend of vegetation changes in this study revealed that the meadow steppe, typical steppe, and alpine steppe exhibited an increasing trend; the desert steppe and the Gobi Desert showed a significant decreasing trend; and only the forest showed no significant trend (Table 1). Tsydypov et al. [22] also found that the variation of the Gobi Desert NDVI decreased in Mongolia and there was a trend of xerophytization had already been occurring in the sparse vegetation of the Gobi Desert. Bao et al. [21] utilized the least-squares linear regression model, and in this study, all pixels of the same vegetation type were included, without considering the spatial heterogeneity of the vegetation changes. We suggest that this is the main reason for the difference in trend results. Global warming has been observed in different studies with different data sources [50][51][52]. According to the records at 48 meteorological stations in Mongolia, Dagvadorj et al. [17] showed that the annual mean temperature of Mongolia increased by 0.062 °C during 1981-2006, which is in agreement with our study. Moreover, Dagvadorj et al. [17] observed that precipitation showed a downward trend, which was also supported by our results. The results of vegetation trend in our findings are different from Zhao et al. [53], as they showed that the NDVI has not shown a persistent increase after 1998, thus causing an insignificant trend over the entire study period from 1982 to 2011. This discrepancy might be caused by the reason that Zhao et al. [53] analyzed the entire Mongolia as a whole to explore the trends of vegetation changes but did not detect the trends by partitioned areas.

Driving Force of Vegetation Distribution and Changes
Several studies have shown that Mongolian vegetation is particularly sensitive and susceptible to climate change and human activity [21][22][23]. Moreover, the spatial distribution of climate factors in Mongolia is mainly because of the combined effects of latitude, topography, and monsoon effects. In this study, we found a significant positive correlation between the spatial distribution of vegetation and the precipitation, which is consistent with the study of Bao et al. [21], Gantsetseg et al. [54], Filei et al. [23], and Zhao et al. [53]. Furthermore, vegetation type and soil type are the main factors affecting vegetation distribution. Note that different vegetation types have different NDVI values. The relationship between soil and vegetation distribution is relatively complicated. From the perspective of soil genesis, the effects of soil on vegetation distribution cannot be simply determined because vegetation also has a transformation effect on soil and different NDVI values (different vegetation types) have different effects on soil transformation [55]. Therefore, there is a relatively strong interaction between soil type and vegetation distribution. Related studies have also shown that when precipitation is the dominant factor, soil type has an important impact on vegetation growth and rainwater reuse [6,56].
Taking the whole of Mongolia as the study area, the influence (q value) of environmental factors on vegetation changes is small. This is because we ignored the spatial heterogeneity of vegetation changes. Therefore, we separately detected the driving force of different vegetation types for their changes during 1982-2015. The results show that grazing was the dominant factor leading to changes in forest, meadow steppe, desert steppe, and Gobi Desert. This indicates that overgrazing had the strongest influence on most of Mongolia's vegetation changes over the studied period. In contrast, for the typical steppe and alpine steppe, topographical factors were dominant because the topography plays an important role in the redistribution of water resources and solar energy resources. For example, the photosynthesis of sunny slope vegetation is stronger than that of shady slopes, and the evapotranspiration of shady slope vegetation is weaker than that of vegetation growing on sunny slopes. Quantifying the driving factors of vegetation distribution has great guiding significance for agricultural site selection. The factor ranking results can provide more accurate factor weight information. In process of constructing the agricultural location model, it can help researchers to assign reasonable weight values to each factor [57].

Optimal Range of Vegetation Growth
Quantification of the optimal range of environmental factors for vegetation distribution in a region helps decision makers to take appropriate interventions to restore vegetation and improve the ecological environment. In GDM, optimal range information is generated from the strata. The division of this strata is a part of the data preprocessing step. When the environmental factor is the continuous numerical data (e.g., precipitation or temperature), the first step is to discretize the environmental factor into strata. The choice of an appropriate discretization method is important because discrete results can directly affect the subsequent GDM results. Many studies recommend using the largest q value as the optimal discretization criterion [58,59]. Therefore, in our study, we used the maximum q value as a criterion function to find the optimal layering method, ensuring the optimal range.
Due to the interaction between the factors, the optimal range of one environmental factor for vegetation growth cannot be considered separately. For example, in current study, the range −12-0 °C was taken as the suitable limit of temperature factor for vegetation growth, however, this range was obtained by combining the temperature with the precipitation and topography factors. The northern part of Mongolia is mainly affected by the airflow of the Arctic Ocean. In addition, the terrain is undulating, and the precipitation is relatively high, causing the dense vegetation. However, because of the high latitude and altitude, the temperature in this area is lower than that of the southern area. The optimal temperature range −12-0 °C for vegetation growth might only reflect the annual average temperature in the northern part of Mongolia.

Relative Degree of Vegetation Changes
The initial value of NDVI can reflect the vegetation coverage at the beginning of vegetation change (i.e., 1982 in this study). Changes in the NDVI are often used as indicators of the degree of desertification development. However, many studies on this aspect only focus on the trend of vegetation change and rarely pay attention to the state of the region before the change [23,60]. To assess the degree of vegetation changes, not only do we need to obtain the trends and the amount of vegetation change, we also need to take the initial state of the vegetation into account. For example, in this study, the NDVI decrease in the Gobi Desert region was between 0 and 0.1 during 1982-2015. However, we should consider that the vegetation index value in this area before the decline was between 0 and 0.2 (initial value).

Caveats of Our Study
Our study is subject to a few caveats. First, we derived the NDVI values from NOAA/AVHRR GIMMS NDVI3g dataset (1/12° and bimonthly frequency, spanning 1982-2015) to explore vegetation distribution and changes. For current study, the spatial resolutions are relatively coarse. As large spatial scale NDVI over heterogeneous surfaces is easily susceptible to soil background in lowvegetated areas and saturated in well-vegetated areas. There are some subsequent vegetation indices to improve it, for example, Enhanced Vegetation Index (EVI) [61] and Soil Adjusted Vegetation Index (SAVI) [62,63]. In addition, biomass estimation can also be used as a good indicator of vegetation conditions. For example, Otgonbayar et al. [64] has estimated pasture biomass in Mongolia using partial least squares, random forest regression, and Landsat 8 imagery. However, considering data availability, this GIMMS NDVI3g dataset is still likely the best one we can currently use for a long time series vegetation analysis. For example, although EVI and SAVI can be derived from a MODIS dataset with higher spatial resolution, MODIS datasets cover a relatively short period of time only since 2000 [65]. Furthermore, long time series vegetation analysis can delve into more objective vegetation change trend information and improve the robustness of analysis results. For vegetation trend research, previous research has shown that the trends of NDVI changes obtained by Landsat (30 m), MODIS (1 km), and AVHRR (8 km) are not significantly different [66].
Second, we used climate data from the Global Land Data Assimilation System (GLDAS) Catchment Land Surface Model (0.25° × 0.25° spatial resolution), which also seem relatively coarse, especially for north of Mongolia (mountainous areas). Obviously, some meteorological data with high temporal-spatial resolution and high quality have been developed in Mongolia, such as average monthly air temperature over Mongolia (spanning 2002-2017) estimated by Otgonbayar et al. [67] using MODIS land surface temperature (LST) time series and machine learning techniques. However, considering the requirement that the meteorological factors have to be resampled to match the spatial resolution of the NDVI data, the GLDAS dataset is a better choice to meet this requirement. In addition, our model results were executed with data at coarse spatial and temporal scales, which possibly result in the potential scale dependence of the model results. However, we are tracing each pixel to calculate the spatial and temporal characteristics of NDVI changes and the corresponding drivers. Additionally, for such study of macro-scale vegetation changes and driving forces, with the improvement of data resolution, the marginal benefits are gradually decreasing. More in-depth research should be further examined on the scale dependence of the trends of NDVI changes.

Conclusions
Based on the GIMMS AVHRR NDVI3g datasets from 1982 to 2015, five climatic data sets (precipitation, air temperature, specific humidity, wind speed, and snow depth) of different time series, four topographical data sets (elevation, curvature, slope degree, and slope aspect), vegetation type data, soil type data, and livestock quantity information, the spatial and temporal distribution characteristic of NDVI changes, the driving force detection of NDVI distribution and changes, the interaction of different environmental factors on vegetation growth, and the optimal range of environmental factors for vegetation growth were analyzed in Mongolia from 1982 to 2015.
On the whole, the area showing a downward trend (accounts for 15% of the total area) was basically equal to the area showing an upward trend (accounts for 14% of the total area) of vegetation variation in Mongolia. However, there was significant heterogeneity between different vegetation types. The meadow steppe, typical steppe, and alpine steppe of Mongolia mainly showed a clear greening tendency, whereas the desert steppe and Gobi Desert showed a degradation tendency.
We found that precipitation, vegetation type, and soil type were the dominant factors influencing NDVI distribution across the whole of Mongolia. The interaction between each pair of factors manifested as mutual enhancement and nonlinear enhancement. Among these environmental factors, topographical factors and climatic factors were most likely to result in nonlinear enhancement effects. This confirms that environmental factors interacted to influence the distribution of NDVI, and interaction between environmental factors plays a more vital role than the influence of the individual one on vegetation growth.
Revealing the optimal range of key environmental factors on vegetation growth is equivalent to clarify the ecological niche of vegetation growth in a study area, facilitating ecological protection and vegetation restoration, and alleviating desertification problems in semi-arid and arid regions. For example, in this study, the results showed that when the annual precipitation is between 331 mm and 596 mm, forest vegetation type and pine sandy soil type could generate the maximum NDVI value in Mongolia.
We stress that the obtained relationships between the environmental factors and vegetation growth are only statistical relationships, rather than causal relationships. However, the results of statistical analysis are often the forerunners of the physical mechanism analysis between geographical elements. Therefore, the results of this study are still instructive to Mongolia's response to climate change in the process of ecological construction and animal husbandry development.

Conflicts of Interest:
The authors declare no conflict of interest.