A Regionalized Study on the Spatial-Temporal Changes of Grassland Cover in the Three-River Headwaters Region from 2000 to 2016

The Three-River Headwaters Region (TRHR) is located in the interior of the Qinghai-Tibetan Plateau, which is a typical research area in East Asia and is of fragile environment. This paper studied the characteristics of grassland cover changes in the TRHR between 2000 and 2016 using methods of area division (AD) based on natural conditions and tabulate area (TA) dependent on Moderate-resolution Imaging Spectroradiometer (MODIS) 44B product. Further investigations were conducted on some of the typical areas to determine the characteristics of the changes and discuss the driving factors behind these changes. Classification and Regression Trees (CART), Random Forest (RF), Bayesian (BAYE), and Support Vector Machine (SVM) Machine Learning (ML) methods were employed to evaluate the correlation between grassland cover changes and corresponding variables. The overall trend for grassland cover in the TRHR towards recovery that rose 0.91% during the 17-year study period. The results showed that: (1) The change in grassland cover was more divisive in similar elevation and temperature conditions when the precipitation was stronger. The higher the temperature was, the more significant the rise of grassland cover was in comparable elevation and precipitation conditions. (2) There was a distinct decline and high change standard deviation of grassland cover in some divided areas, and strong correlations were found between grassland cover change and aspect, slope, or elevation in these areas. (3) The study methods of AD and TA achieved enhancing performance in interpretation of grassland cover changes in the broad and high elevation variation areas. (4) RF and CART methods showed higher stability and accuracy in application of grassland cover change study in TRHR among the four ML methods utilized in this study.


Introduction
The Three-River Headwaters Region (TRHR) is a typical study area with high ecology values [1][2][3].The main vegetation cover type in TRHR is grassland, which conserves water resource, preserves soil, and grazes livestock [4][5][6].The grassland cover is traditionally studied by in-situ measurement, but it is inefficient in grassland cover change monitoring due to a large area of grassland [7,8].Thus, the remote sensing images, product indicators, and methods were much more effective in researches of grassland cover change [9,10].Schmidt et al. [11] calculated cover and management factor (C-factor) by a combination of high spatial resolution Swissimage false-color infrared (Swissimage FCIR) and high temporal resolution Fraction of green Vegetation Cover (FCover300m) in Switzerland.Wang et al. [12] examined the correspondence between Soil Moisture Content (SMC) and vegetation cover via MODIS13A3 Normalized Difference Vegetation Index (NDVI) product on Loess Plateau, while Liu et al. [13] studied the vegetation variation through Global Inventory Modeling and Mapping Studies (GIMMS) NDVI 3g product and Theil-Sen median analysis in Xinjiang.Li et al. [14] monitored and modeled the grassland dynamics by MODIS13Q1 product and Mann-Kendall trend analysis over accumulated precipitation and stocking intensity factors in Saskatchewan, Canada.Li et al. [15] monitored and analyzed the grassland desertification using visual interpretation results of Landsat TM/ETM+ images of 1993, 2000, 2006, and 2011 in Ningxia Hui Autonomous Region, China.Si et al. [16] validated the quantity and quality of grassland through Leaf Area Index (LAI) and Canopy Chlorophyll Content (CCC) indicators from Medium Resolution Imaging Spectrometer (MERIS) over Groningen and Friesland provinces in Netherland.
Meanwhile, there has been a large number of previous studies of grassland cover changes on the TRHR.Liu et al., Shao et al., and Liu et al. [4][5][6] studied the spatial-temporal changes in TRHR via comparison of visual interpretation results of Landsat MSS image data in late 1970's, TM data in early 1990's, TM/ETM data in 2004, and Huan Jing (means "environment", HJ) satellite data in 2012.The spatial resolution of the interpretation results was high, while it could be a time and effort consuming work.Liu et al. [1] analyzed spatial-temporal characteristics of MODIS13Q1 NDVI on county and basin scales between 2000 and 2011, and they derived the changing tendencies of grassland cover by Hurst Index without considering the relationship between the tendencies and terrain.Shen et al. [2] studied vegetation changes in TRHR through MODIS09GA daily NDVI data and Detecting Breakpoints and Estimating Segments in Trend (DBEST).However, they did not take the spatial heterogeneity of natural conditions into consideration in TRHR, which might lead to inaccuracy of the DBEST analysis.
To sum up, the previous studies of grassland cover established the models or processed the analysis across the entire study areas.However, there would be observable spatial heterogeneity on natural conditions when the study area was vast, which might result in the accuracy loss of the results.As TRHR is a broad area with high elevation range and complicate environment, the area division (AD) based on natural conditions before analysis process may reduce the spatial heterogeneity and achieve better results.The AD also fit the ecological areas classification principals of climate and geomorphology differences within different areas [17].Meanwhile, commonly the previous studies researched the changes of grassland cover over time series (annually, seasonally, monthly, or daily), while the characteristics of relationship between grassland cover changes and terrain (elevation (DEM), slope (SLP), or aspect (ASP)) were seldom extended.The TA method applied in this study revealed the relationship between grassland cover changes and terrain factors, and the ML algorithms of CART, RF, BAYE, and SVM provided another way to indicate the properties of the relationship.
The objectives of this study were to derive the characteristics of relationship between grassland cover changes and terrain factors (DEM, SLP, or ASP), and evaluate the characteristics through CART, RF, BAYE, and SVM machine learning methods.Natural condition factors of DEM, air temperature (TEM), and precipitation (PRE) data were applied to divide the TRHR, so the natural conditions could be more homogeneity in each divided areas.The Percent Non Tree Vegetation (NTV) band of MODIS44B product was used as the grassland cover indicator.The Tabulate Area (TA) method was brought in to derive the relationship between NTV changes and terrain factors (DEM, SLP, or ASP).Moreover, the CART, RF, BAYE, and SVM ML algorithms were utilized to evaluate the relationship between terrain factors and the characteristics derived from TA method.

Study Area
The TRHR is situated in the Northeast of the Qinghai-Tibetan Plateau (QTP) and is the source of the Yangtze, Yellow, and Lancang Rivers.The area of the TRHR covers more than 3.3 × 10 5 km 2 and is commented as the "Chinese Water Tower" [4].The TRHR is bounded by 31 • N-37 • N latitude and 89 • E-103 • E longitude.The elevation of the west part is higher than the east, which makes the west part dryer and colder, as shown in Figure 1.The climate of TRHR is arid and frigid owe to high elevation [1][2][3][4][5][6], and it has changed in the TRHR due to global warming processes [17][18][19][20][21][22][23][24].The TEM in the QTP has increased 0.6 • C/10 years in the past 30 years, which is almost twice the rate from the past 50 years.Meanwhile, the PRE in the same area has risen 1.4 mm/10 years during the last 30 years, the PRE ascending velocity is almost one sixth of the rate during the past 50 years [18].Thus, there are significant warm-drying trends in the warm-wetting processes.One observation showed that the land surface temperature (LST) in the QTP rose 1.8 • C during the past 50 years [19].Simultaneously, human activities such as the raising livestock, construction, and the Conservation and Restoration Project of grassland in TRHR led to the turbulence in the changes over the entire region [20], these changes established a deep influence in the TRHR, QTP, and even East Asia [21][22][23][24][25].

Data Source and Preprocessing
The NTV data is a layer of the MOD44B Vegetation Continuous Fields product which is composited by Terra MODIS 250 m and 500 m land surface reflectance data.MOD44B product provides 250 m yearly global retrieval vegetation cover data from 2000.The radiation and cloud cover influences were eliminated to ensure the reliability of product data [26,27].Different from the previous vegetation cover products, the product contained the NTV, percent tree cover, and percent non vegetated bands.The values ranged from 1% to 100%, and pixels were identified as water when the value equaled 200.In this study, we selected 17 years of NTV data from 2000 to 2016 and removed the water pixels, the TRHR, and divided areas (in Section 4.1) border mask layers were applied to extract the study areas.
Gessner et al. [28] made a validation of the MOD44B product through higher spatial resolution remote sensing data (QuickBird, 0.6 m-1 m and Landsat, 30 m) and in-situ measurement; the results showed that the Root Mean Square Deviation (RMSE) varied from 3.1% to 8.2% when compared with higher resolution data, while the RMSE ranged between 4.4% and 9.9% when compared with in-situ measurement.Therefore, they considered the high accuracy of the MOD44B product and the high value of its applications.
In this study, we utilized the Chinese monthly surface precipitation 0.5 • × 0.5 • grid dataset from 2000 to 2015 and the Chinese monthly surface temperature 0.5 • × 0.5 • grid dataset from 2000 to 2014, which were provided by the National Atmospheric Science Data Sharing and Serving Infrastructure (http://data.cma.cn/site/index.html).The PRE and TEM data were produced by the Infrastructure through ANUSPLIN from meteorological station data [29,30], together with Shuttle Radar Topography Mission (SRTM) DEM data [31] with a spatial resolution of 90 m, the resampling of the PRE, TEM, and DEM raster data through cubic convolution interpolation was processed to ensure the spatial resolution was equal to the NTV data.By calculating the DEM data, we were able to achieve the slope and aspect data needed for the TA method of TRHR.

Area Division Process
In this study, we attempted to divide the TRHR based on natural conditions.The entire region was divided into research areas with similar natural conditions of DEM, TEM, and PRE.We then

Area Division Process
In this study, we attempted to divide the TRHR based on natural conditions.The entire region was divided into research areas with similar natural conditions of DEM, TEM, and PRE.We then investigated the NTV variations (NTVV) and correlate drivers in areas with similar natural conditions.
The mean PRE raster data over 16 years and mean TEM data over 15 years were calculated by resampled ANUSPLIN monthly data.Approximate bisection of 476.88 mm for PRE and −3.31 • C for TEM values were used to divide the precipitation and temperature data, while section points 4400 m and 4900 m were utilized to divide the DEM data, as shown in Table 1.The union overlay analysis was processed over the three data layers described in Table 1, and resulted in the 12 divided areas shown in Figure 2. The mean PRE raster data over 16 years and mean TEM data over 15 years were calculated by resampled ANUSPLIN monthly data.Approximate bisection of 476.88 mm for PRE and −3.31 °C for TEM values were used to divide the precipitation and temperature data, while section points 4400 m and 4900 m were utilized to divide the DEM data, as shown in Table 1.The union overlay analysis was processed over the three data layers described in Table 1, and resulted in the 12 divided areas shown in Figure 2. The research area classifications are represented by a three-letter combination code: The first letter classifies the elevation (A, B, C); the second letter is the temperature category (A, B); and the third letter represents the precipitation class (A, B).For example, AAA represent a research area where the elevation is lower than 4400 m, the temperature lower than -3.31 °C , and the precipitation is lower than 476.88 mm, as further described in Table 2.The research area classifications are represented by a three-letter combination code: The first letter classifies the elevation (A, B, C); the second letter is the temperature category (A, B); and the third letter represents the precipitation class (A, B).For example, AAA represent a research area where the elevation is lower than 4400 m, the temperature lower than −3.31 • C, and the precipitation is lower than 476.88 mm, as further described in Table 2.In following sections, we selected the most significantly decreased area AAB and the areas with higher standard variations (BAB, BBB, CAB, and CBB) for further analysis.Among those areas, BAB and BBB covered large regions in the east and west directions.Due to the uncertainty brought in by the complicated natural conditions and varied climate conditions in the plateau, we divided the two areas into East, Center, and West parts, marked as (I), (II), and (III).

Tabulate Area Method
The TA algorithm is used to count some statistics (mean value, median value, standard deviation, et al.) in two raster images with same resolution.One of the images is considered as the zonal data to divide the statistics area in the other image, while the other image is the statistic data which implement statistics calculations based on areas determined by prior image [32].In this study, the mean value was selected to demonstrate the variation rules of grassland cover change, as shown in Figure 3. Firstly, the pixels with same value in data Figure 3(b) are determined, and then the same pixels located in Figure 3(a) are used to calculate the mean value of Figure 3(a) to the establish raster shown in Figure 3(c).
In this research, we processed NTVV data though TA algorithm by DEM, SLP and ASP data.The slope and aspect data values were transformed to integer according to the rules of TA algorithm.

Machine Learning Methods
We utilized CART, RF, BAYE, and SVM ML regression algorithms to evaluate the NTVV results based on study AD and TA methods.Brief introductions of the algorithms listed as follows.
CART was a top-down classify process which divided the whole dataset into two parts, a binary decision tree [33].The division process resulted in the maximum homogeneity in the same divided dataset.Then the process continued until it reached the threshold set by the user or the dataset could not be divided any more.The dataset would be called leave node then.The prediction was made based on the nodes [34].
RF was established by Leo Breiman and Adele Cutler to achieve better prediction results through multiple decision trees [35].A random forest classifier was set up, and several sub-sets were divided from the data first.Then the amount of sub-sets increased to a number defined by the user, each one of the origin data was classified into a certain tree.The prediction could be made according to the individual trees [36].
There were two main steps in BAYE algorithm: Firstly, the likelihood function of training dataset was obtained, we calculated the posterior distribution through the relationship of prior probability and the likelihood function.Then, the posterior distribution was utilized to make classification of the dataset to be predicted [37].
SVM was improved by Vapnik and it was widely used in the ML model in text categorization, image classification and many other fields [38].It was used to classify data into one of two or more categories.The major principle of the SVM was to keep the largest distance between two categories.When the dataset was not linearly separable, it could be mapped into higher-dimensional space.
A 500 m × 500 m dot matrix was utilized to sample the NTVV, DEM, ASP, and SLP data in further studied areas.The samples were used as the training set to estimate the NTVV of original scale 250 m, the amount of training set was about 25%.We applied RF, CART, BAYE, and SVM algorithms to make prediction of NTVV based on corresponding variables in individual areas, ASP for AAB, DEM for BAB(I), BAB(II), BBB(I), BBB(II), and SLP for BAB(III), BBB(III), CAB, and CBB.Then, we made comparison between the predicted data and the MOD44B NTVV data processed by the TA algorithm.Meanwhile, the 500 m × 500 m dot matrix training set was used to estimate NTVV in whole TRHR based on individual and all three independent variables without AD and TA algorithms.The comparison between ML results and the AD and TA methods results was made, in order to evaluate the efficiency of AD and TA methods applied in this study.Moreover, the validation of TA method applied in divided areas was made through ML methods.The package scikit-learn was Firstly, the pixels with same value in data Figure 3b are determined, and then the same pixels located in Figure 3a are used to calculate the mean value of Figure 3a to the establish raster shown in Figure 3c.
In this research, we processed NTVV data though TA algorithm by DEM, SLP and ASP data.The slope and aspect data values were transformed to integer according to the rules of TA algorithm.

Machine Learning Methods
We utilized CART, RF, BAYE, and SVM ML regression algorithms to evaluate the NTVV results based on study AD and TA methods.Brief introductions of the algorithms listed as follows.
CART was a top-down classify process which divided the whole dataset into two parts, a binary decision tree [33].The division process resulted in the maximum homogeneity in the same divided dataset.Then the process continued until it reached the threshold set by the user or the dataset could not be divided any more.The dataset would be called leave node then.The prediction was made based on the nodes [34].
RF was established by Leo Breiman and Adele Cutler to achieve better prediction results through multiple decision trees [35].A random forest classifier was set up, and several sub-sets were divided from the data first.Then the amount of sub-sets increased to a number defined by the user, each one of the origin data was classified into a certain tree.The prediction could be made according to the individual trees [36].
There were two main steps in BAYE algorithm: Firstly, the likelihood function of training dataset was obtained, we calculated the posterior distribution through the relationship of prior probability and the likelihood function.Then, the posterior distribution was utilized to make classification of the dataset to be predicted [37].
SVM was improved by Vapnik and it was widely used in the ML model in text categorization, image classification and many other fields [38].It was used to classify data into one of two or more categories.The major principle of the SVM was to keep the largest distance between two categories.When the dataset was not linearly separable, it could be mapped into higher-dimensional space.
A 500 m × 500 m dot matrix was utilized to sample the NTVV, DEM, ASP, and SLP data in further studied areas.The samples were used as the training set to estimate the NTVV of original scale 250 m, the amount of training set was about 25%.We applied RF, CART, BAYE, and SVM algorithms to make prediction of NTVV based on corresponding variables in individual areas, ASP for AAB, DEM for BAB(I), BAB(II), BBB(I), BBB(II), and SLP for BAB(III), BBB(III), CAB, and CBB.Then, we made comparison between the predicted data and the MOD44B NTVV data processed by the TA algorithm.Meanwhile, the 500 m × 500 m dot matrix training set was used to estimate NTVV in whole TRHR based on individual and all three independent variables without AD and TA algorithms.The comparison between ML results and the AD and TA methods results was made, in order to evaluate the efficiency of AD and TA methods applied in this study.Moreover, the validation of TA method applied in divided areas was made through ML methods.The package scikit-learn was employed and proper parameters were chosen to achieve better results [39].R 2 , RMSE and Mean Absolute Error (MAE) were used to assess the ML results.The specific study steps were shown in Figure 4. employed and proper parameters were chosen to achieve better results [39].R 2 , RMSE and Mean Absolute Error (MAE) were used to assess the ML results.The specific study steps were shown in Figure 4.

Overall Characteristics of the TRHR
In Figure 5, the NTVV between 2000 and 2016 in the TRHR was shown.In 17 years, the NTV increased 0.91%, from 47.26% to 48.17%.We classified the NTVV of TRHR in 17 years into five categories: The serious recession area (NTVV < −15%) was 22,107 km 2 , the mild recession area (−15% ≤ NTVV < −5%) was 53,191 km 2 , the insignificantly changed area (−5% ≤ NTVV < 5%) was 173,080 km 2 , the slightly increased area (5% ≤ NTVV < 15%) was 57,333 km 2 , and the significantly risen area (NTVV ≥ 15%) was 27,401 km 2 .By calculating the annual average PRE raster data for 16 years, the TEM data over 15 years in the TRHR, and the NTV data over 17 years for the linear regression analysis, the results showed that the R 2 value of the NTV was 0.144.The upward tendency was 0.98%/10a; a minor increasing rate compared with the research of Shao et al. [40].The R 2 of the PRE was 0.09, while the increasing rate was 35.8 mm/10a; values close to those obtained by Li et al. [41].The R 2 of the TEM was 0.4, while the upward trend was 0.66 °C /10a; a milder result than obtained by Yi et al. [42].
The results of AD show that the NTV decreased in both low-elevation and low-temperature areas (AAA, AAB), both high-elevation and low-temperature areas (CAA, CAB), and the highelevation high-temperature and high-precipitation area (CBB).The NTV increased in both lowelevation and high-temperature areas (ABA, ABB), the moderate-elevation area and high-elevation high-temperature and low-precipitation area (CBA).Among these areas, the NTV in the AAA, AAB, BBB, and CBA areas increased more significantly.In areas with similar elevation and temperature, the more the precipitation was, the higher the standard deviation was, meaning that the NTV changed in a more complicated way than anticipated.

Area with Significantly NTV Decrease
In the AAB area, the decrease phenomenon concentrated in the western part of the Amne Machin Mountain, southwest of Xinghai County, northwest of Maqeen County.Together with the aspect data, the results showed that pixels in the West (247.5°-292.5°),Northwest (292.5°-337.5°),North (337.5°-22.5°),and Northeast (22.5°-67.5°)slopes, representing 49.19% of all the pixels, had NTV that decreased 3.51% in 17 years, while the NTV in other aspects recessed 0.41%, as shown in Figure 6b.By calculating the annual average PRE raster data for 16 years, the TEM data over 15 years in the TRHR, and the NTV data over 17 years for the linear regression analysis, the results showed that the R 2 value of the NTV was 0.144.The upward tendency was 0.98%/10a; a minor increasing rate compared with the research of Shao et al. [40].The R 2 of the PRE was 0.09, while the increasing rate was 35.8 mm/10a; values close to those obtained by Li et al. [41].The R 2 of the TEM was 0.4, while the upward trend was 0.66 • C/10a; a milder result than obtained by Yi et al. [42].
The results of AD show that the NTV decreased in both low-elevation and low-temperature areas (AAA, AAB), both high-elevation and low-temperature areas (CAA, CAB), and the high-elevation high-temperature and high-precipitation area (CBB).The NTV increased in both low-elevation and high-temperature areas (ABA, ABB), the moderate-elevation area and high-elevation high-temperature and low-precipitation area (CBA).Among these areas, the NTV in the AAA, AAB, BBB, and CBA areas increased more significantly.In areas with similar elevation and temperature, the more the precipitation was, the higher the standard deviation was, meaning that the NTV changed in a more complicated way than anticipated.

Area with Significantly NTV Decrease
In the AAB area, the decrease phenomenon concentrated in the western part of the Amne Machin Mountain, southwest of Xinghai County, northwest of Maqeen County.Together with the aspect data, the results showed that pixels in the West (247.5 • -292.5 • ), Northwest (292.5 • -337.5 • ), North (337.5 • -22.5 • ), and Northeast (22.5 • -67.5 • ) slopes, representing 49.19% of all the pixels, had NTV that decreased 3.51% in 17 years, while the NTV in other aspects recessed 0.41%, as shown in Figure 6b.There was an observable correlation between the NTVV and aspect in area AAB.The NTV decreased significantly on the NE, N, NW, and W slopes, while the NTV on other aspects remained stable.The area was a populous place in the Guolo autonomous prefecture with a strong human activity influence and livestock, together with rapidly rising precipitation, resulting in the recession of the NTV.The more notable regression on the ubac slope might ascribe to the thinner grass density, more fragile ecology that is more sensitive to natural condition changes and weaker resistibility to wind and rain erosion [25,[43][44][45][46][47].

Areas with High NTVV Standard Deviation
Compared with the NTV significantly decreased area AAB, the elevation of the complex variation areas was higher and the natural conditions more complicated, which led to a more varied NTVV.
Subarea I of area BAB was located at center of the Guoluo autonomous prefecture and its surrounding areas lay between the Amne Machin Mountain and the Bayan Har Mountain, as shown in Figure 7a.The NTVV was −6.43% during the 17-year study period and the recession phenomenon was significant.Processing of the TA along with the elevation data showed that the NTV rose in areas higher than 4800 m, while decreased at lower elevations.The most significantly decreased area was located at about 4500 m and achieved −16~−17%, as shown in Figure 7b.Subarea II was situated at the center of the TRHR, bordering Madoi and Chindu counties, as shown in Figure 7c.The NTVV was 1.86% during the 17-year study period which was inversely proportional to the elevation between 4450 m and 4830 m.The NTV increased between 4430 m and 4720 m and recessed over 4720 m, as illustrated in Figure 7d.Subarea III was located in the East of Zhidoi County, Tongtian River Basin, as shown Figure 7e.The NTV observably increased by 4.97% in the 17-year study period, and through the processing of the TA algorithm to slope, the results showed that the NTV insignificantly changed between 20° and 25°, while it increased in areas with other slopes, as shown in Figure 7f.There was an observable correlation between the NTVV and aspect in area AAB.The NTV decreased significantly on the NE, N, NW, and W slopes, while the NTV on other aspects remained stable.The area was a populous place in the Guolo autonomous prefecture with a strong human activity influence and livestock, together with rapidly rising precipitation, resulting in the recession of the NTV.The more notable regression on the ubac slope might ascribe to the thinner grass density, more fragile ecology that is more sensitive to natural condition changes and weaker resistibility to wind and rain erosion [25,[43][44][45][46][47].

Areas with High NTVV Standard Deviation
Compared with the NTV significantly decreased area AAB, the elevation of the complex variation areas was higher and the natural conditions more complicated, which led to a more varied NTVV.
Subarea I of area BAB was located at center of the Guoluo autonomous prefecture and its surrounding areas lay between the Amne Machin Mountain and the Bayan Har Mountain, as shown in Figure 7a.The NTVV was −6.43% during the 17-year study period and the recession phenomenon was significant.Processing of the TA along with the elevation data showed that the NTV rose in areas higher than 4800 m, while decreased at lower elevations.The most significantly decreased area was located at about 4500 m and achieved −16~−17%, as shown in Figure 7b.Subarea II was situated at the center of the TRHR, bordering Madoi and Chindu counties, as shown in Figure 7c.The NTVV was 1.86% during the 17-year study period which was inversely proportional to the elevation between 4450 m and 4830 m.The NTV increased between 4430 m and 4720 m and recessed over 4720 m, as illustrated in Figure 7d.Subarea III was located in the East of Zhidoi County, Tongtian River Basin, as shown Figure 7e.The NTV observably increased by 4.97% in the 17-year study period, and through the processing of the TA algorithm to slope, the results showed that the NTV insignificantly changed between 20 • and 25 • , while it increased in areas with other slopes, as shown in Figure 7f.Further study of subarea I in area BAB indicated that the NTV decreased most near to 4500 m.Considering the main settlement (Changmahe Township, 4460 m) nearby, we might argue that the human activities were intensive near that elevation.Human activities decreased with the rise of Further study of subarea I in area BAB indicated that the NTV decreased most near to 4500 m.Considering the main settlement (Changmahe Township, 4460 m) nearby, we might argue that the human activities were intensive near that elevation.Human activities decreased with the rise of elevation.Thus, the NTVV rose linearly with elevation.The increase of NTV above 4800 m might ascribe to rare human activities and a warm-wetting tendency above 4800 m.In subarea II, the relationship between the NTVV and elevation was quite different from subarea I where the NTV rose beneath 4700 m.Considering that the subarea was located in the middle region of the Bayan Har Mountain where the Bayan Har Mountain pass was found.There were few human activities, together with the warm-wetting tendencies in the TRHR, which might result in the rise of the NTV.The Western subarea III was located in the midfields of the Bayan Har and Tanggula Mountains where the parallel mountain gathers the warm-wetting clouds, which might result in the increase of the NTV.The erosion led to the insignificant changes between 20 • and 25 • in subarea III [48][49][50].
The NTV observably rose during the 17-year study period in area BBB, where the standard variation was the highest among the divided study areas.Subarea I was distributed in Maqeen, Gadee, and Tarlag Counties, the middle and east parts of the Amne Machin Mountain, the center of Bayan Har Mountain, as shown in Figure 8a In subarea I of BBB the curve of the NTVV with elevation was similar to subarea BAB(I), but the value field was higher.The curve was similar to subarea BAB(II) in sub area II, but was smoother.The areas with higher temperature had a more stable NTVV in similar precipitation and elevation areas.Subarea III, which was located at east of Tanggula Mountain, had better hydrothermal conditions due to the mountain blocked the warm and wet flow.With the rose of the slope, the erosion influenced the ascent of the NTV and the Conservation and Restoration Project activities might positively affect the grassland cover changes as well [25].
The NTV in area CAB remained almost steady during the 17-year study period.The area was distributed in the West of Zadoi County, Southwest of the Tanggula Township of Golmud City, as shown Figure 9a.The NTV decreased below 11°, while the NTV increased above 11°, as demonstrated in Figure 9b.The NTVV rose with the slope, and the NTV decreased in low slope areas, it might ascribe to the permafrost thaw caused by rises in the air temperature [51][52][53], the grass erosion caused by pests and In subarea I of BBB the curve of the NTVV with elevation was similar to subarea BAB(I), but the value field was higher.The curve was similar to subarea BAB(II) in sub area II, but was smoother.The areas with higher temperature had a more stable NTVV in similar precipitation and elevation areas.Subarea III, which was located at east of Tanggula Mountain, had better hydrothermal conditions due to the mountain blocked the warm and wet flow.With the rose of the slope, the erosion influenced the ascent of the NTV and the Conservation and Restoration Project activities might positively affect the grassland cover changes as well [25].
The NTV in area CAB remained almost steady during the 17-year study period.The area was distributed in the West of Zadoi County, Southwest of the Tanggula Township of Golmud City, as shown Figure 9a.The NTV decreased below 11 • , while the NTV increased above 11 • , as demonstrated in Figure 9b.In subarea I of BBB the curve of the NTVV with elevation was similar to subarea BAB(I), but the value field was higher.The curve was similar to subarea BAB(II) in sub area II, but was smoother.The areas with higher temperature had a more stable NTVV in similar precipitation and elevation areas.Subarea III, which was located at east of Tanggula Mountain, had better hydrothermal conditions due to the mountain blocked the warm and wet flow.With the rose of the slope, the erosion influenced the ascent of the NTV and the Conservation and Restoration Project activities might positively affect the grassland cover changes as well [25].
The NTV in area CAB remained almost steady during the 17-year study period.The area was distributed in the West of Zadoi County, Southwest of the Tanggula Township of Golmud City, as shown Figure 9a.The NTV decreased below 11°, while the NTV increased above 11°, as demonstrated in Figure 9b.The NTVV rose with the slope, and the NTV decreased in low slope areas, it might ascribe to the permafrost thaw caused by rises in the air temperature [51][52][53], the grass erosion caused by pests and The NTVV rose with the slope, and the NTV decreased in low slope areas, it might ascribe to the permafrost thaw caused by rises in the air temperature [51][52][53], the grass erosion caused by pests and rodents [54], and the "black soil beach" affects [55,56].The fall of temperature caused by increase of elevation, and the melt water and weaken pest effects caused by farther distance from water with the improvement of the slope angle, together with the overall warm-wetting tendency, might be the reason of the increase of the NTV in area CAB.
The NTV recessed lightly during the 17-year study period in the CBB area, which is mainly situated at the Center East of Zadoi County and the Eastern part of the Tanggula Mountain, as shown in Figure 10a.The NTV decreased below 19 • , while it rose above this slope, as illustrated in Figure 10b.
Sustainability 2018, 10, x FOR PEER REVIEW 14 of 23 rodents [54], and the "black soil beach" affects [55,56].The fall of temperature caused by increase of elevation, and the melt water and weaken pest effects caused by farther distance from water with the improvement of the slope angle, together with the overall warm-wetting tendency, might be the reason of the increase of the NTV in area CAB.
The NTV recessed lightly during the 17-year study period in the CBB area, which is mainly situated at the Center East of Zadoi County and the Eastern part of the Tanggula Mountain, as shown in Figure 10a.The NTV decreased below 19°, while it rose above this slope, as illustrated in Figure 10b.Area CBB was located at the northern slope in the Eastern Tanggula Mountain.The area was warmer than CAB, the warm-drying phenomenon was more significant.Thus, the NTVV value in CBB was lower than CAB.In the high elevation areas CAB and CBB, compared with lower areas, the differences of temperature rise were insignificant, while the differences of precipitation rise were more remarkable during the study period, and this might be the one of factors for the NTVV divergences among high elevation areas.

Machine Learning Results Evaluation
To evaluate the AD and TA methods applied in this study, and to validate the performance of four ML methods selected, the results of comparison between predicted NTVV by ML algorithms and the TA results of NTVV were demonstrated in this section.There were three types of terrain factors, ASP for AAB, DEM for BAB(I), BAB(II), BBB(I), BBB(II), and SLP for BAB(III), BBB(III), CAB, CBB.All the results and corresponding scatterplots, linear regression equations, R2, RMSEs, and MAEs were showed in Figure 11.Area CBB was located at the northern slope in the Eastern Tanggula Mountain.The area was warmer than CAB, the warm-drying phenomenon was more significant.Thus, the NTVV value in CBB was lower than CAB.In the high elevation areas CAB and CBB, compared with lower areas, the differences of temperature rise were insignificant, while the differences of precipitation rise were more remarkable during the study period, and this might be the one of factors for the NTVV divergences among high elevation areas.

Machine Learning Results Evaluation
To evaluate the AD and TA methods applied in this study, and to validate the performance of four ML methods selected, the results of comparison between predicted NTVV by ML algorithms and the TA results of NTVV were demonstrated in this section.There were three types of terrain factors, ASP for AAB, DEM for BAB(I), BAB(II), BBB(I), BBB(II), and SLP for BAB(III), BBB(III), CAB, CBB.All the results and corresponding scatterplots, linear regression equations, R2, RMSEs, and MAEs were showed in Figure 11.The predicted NTVV in area AAB based on ASP showed the lowest R 2 values, while the RMSE and MAE values were at the medium level among the areas.The estimated NTVV in area BAB(I), BAB(II), BBB(I) and BBB(II) by DEM showed the R 2 value ranged 0.032-0.885,while the range of RMSE was 1.846-4.672%,MAE ranged from 1.240% to 4.128%.The NTVV predicted by SLP in area BAB(III), BBB(III), CAB, and CBB illustrated R 2 value between 0.225 and 0.947, the value of RMSE and MAE of RF and CART were below 1%, the details were shown in Table 4.It could be realized the predictions made by slope was more accurate than ASP and DEM, this could be explained by the concentration of the sample.There were dozens of variable values of SLP, while the amount of variable values of DEM and ASP was counted by hundred.Thus, there were more samples for each of the corresponding dependent variable of SLP than elevation and aspect, which made the result more accurate.

Evaluation of Machine Learning Results without AD and TA Methods
In order to validate the effect of AD and TA methods applied in this study, we estimated 250 m NTVV through DEM, SLP, and ASP variables from the 500 m × 500 m dot matrix training set without AD and TA progresses.Then compared it with the ML analysis results in Section 4.4 to evaluate the methods applied in this study, the details of the ML results with AD and TA methods were shown in Table 4.The results of four ML methods without AD and TA in TRHR by DEM, SLP, ASP, or all three terrain factors showed that the R 2 values were less than 0.005, RMSE values were higher than 20%, and MAE values were all beyond than 13%.The results were far less accurate than the ML results with AD and TA methods.
Across the four ML algorithms applied in this study, the SVM method always led to better R 2 values, but the RMSE and MAE values were unstable, there was even unexpected low R 2 value in area BBB(III).The R 2 of BAYE were frequently the lowest among the methods, and unexpected low values occurred in AAB and BBB(I).The RMSE and MAE values of BAYE were sometimes lower than the other methods.Comparing with SVM and BAYE methods, the performance of RF and CART were not always better in the prediction of NTVV.However, there was no abnormal R 2 , RMSE, or MAE value through RF and CART methods.Moreover, the RF and CART methods always performed better in prediction of NTVV on slope.Therefore, RF and CART methods might be more appropriate in application of grassland cover change studies in TRHR with higher stability and accuracy.
It was observably showed that the R 2 , RMSE and MAE results of further studied areas were much better than counterpart of TRHR without AD and TA progresses.The ML results of divided areas based on single variable achieved superior performance than the result of TRHR based on single or all three variables.Thus, the AD and TA study methods used in this study was effective in the interpretation of NTVV.

Evaluation of Machine Learning Results without TA Method
The AD and TA methods improved the interpretation of NTVV over terrain factors, but whether the TA method worked in the divided areas was still a question.The 250 m NTVV was estimated in further studied divided areas by the best performed RF ML method through all three terrain factors (ASP, DEM, and SLP).Together with the RF ML results in divided areas of Sections 4.4 and 4.5, the validation between estimated NTVV and TA method results was made, as shown in Table 5.The R 2 values of RF ranged 0.185-0.949with TA method in further studied areas, while the values ranged from 0.000 to 0.075 without TA.The RMSE ranged 0.593-4.337with TA method, while ranged from 23.092 to 30.145 without TA.The MAE values ranged 0.399-3.388and 17.323-22.676with and without TA separately.It could be derived that the TA method remarkably improved the R 2 , RMSE, and MAE indicators performance of RF ML results.It means that the TA method results were much closer to the RF ML results than the NTVV without TA method.Therefore, the TA method might help for the improvement of NTVV interpretation over single terrain factor in divided areas.

Discussion
The overall NTV increased 0.91% in the TRHR during the 17-year study period, but the NTVV in the different study areas varied.The Conservation and Restoration Project succeed in these areas, such as low elevation northeast part of TRHR, area between north of Bayan Har Mountain and Gyaring and Ngoring Lake, and the region between Tanggula and Bayan Har Mountain, except for the low elevation Northeast part of TRHR, the other areas located in the project area.In area BAB(III), the NTV increased 4.97% during the 17-year study period.However, the NTV decreased at Amne Machin Mountain, Bayan Har Mountain, and the East Tanggula Mountain.The higher the DEM in these areas is, the more vulnerable the environment is, and the NTV deceased 6.43% in area BAB(I) during the 17-year study period.
In areas with similar DEM and PRE, the NTVV showed more positive trends when the TEM was higher, while in areas with similar DEM and TEM, the NTVV showed more variable tendencies with the rise of PRE.
It was observably illustrated that the NTVV had the tendency to follow trends related to some of the factors after the study AD based on natural conditions.The NTV demonstrated distinct correlations to the natural aspects in area AAB, and the NTV in the adret showed insignificant changes while it was evidently recessed in other aspects.
Simultaneously, in NTV complicated changed areas, the NTVV showed inversely proportional responses to elevations below 4500 m while directly proportional responses were observed above in area BAB(I).The NTVV demonstrated adversely proportional responses to elevations between 4450 m and 4830 m in area BAB(II).In area BAB(III), the NTVV decreased with the rise in slope beneath 22 • and increased over 22 • .The NTVV in area BBB(I) was inversely proportional to elevations below 4530 m, while directly proportional above 4530 m.In area BBB(II), the NTVV rose with the rise of elevations beneath 4750 m, and decreased with elevation increases above 4800 m.The NTVV in area BBB(III) had an adversely proportional responses to slope, while it showed directly proportional responses to the slopes in areas CAB and CBB.
In NTV complicated changed areas, the NTV decreased 6.43% during the 17-year study period in area BAB(I), recessed below 4800 m, and most significantly at about 4500 m.In area BAB(II), the NTV increased 1.86% and showed an upward trend between 4430 m and 4720 m, while it decreased above 4720 m.The NTV in area BAB(III) improved 4.97%, especially in slopes between 0 • and 20 • , and also above 25 • .The NTV in area BBB(I) recessed 3.51% and showed a decreasing tendency between elevations of 4400-4730 m and decreased the most at about 4530 m.The NTV increased between 4400 m and 4640 m, while it recessed over 4640 m in area BBB(II).The NTV significantly increased 8.31% in area BBB(III).The NTV in area CAB slightly decreased by 0.33%, it recessed below 11 • , and rose above.In area CBB, the NTV decreased 0.74%, recessed beneath 19 • , and enhanced above 19 • .
The comparison between ML results in divided areas through TA and ML results in TRHR without TA indicated the study methods of AD and TA utilized were efficient in interpretation of grassland cover changes in TRHR.
There were studies on grassland cover on TRHR through visual interpretation results of Landsat and HJ satellite data [4][5][6].The land cover types and land cover change were properly analyzed.The percentages of land cover and cover changes were demonstrated.If the TA method was utilized, the relationship between land cover changes and relative factors (PRE, TEM, DEM, et al.) could be further derived.In a broad region like TRHR, where the area is more than 3.3 × 10 5 km 2 with high DEM value range and complex environment.The accuracy of overall trending analysis [1,2] might not always achieve the application demands.If the AD by related factors (DEM, land cover, climate, et al.) method was employed, the corresponding thresholds and factors value range in the analysis process could be limited, which might lead to better results.The grassland cover change could be more accurately described by TA method instead of depicting by the position of the county or basin.
The grassland cover change was a complicated process and it was difficult to be explained by a single independent variable.In this study, an AD method was proposed to make the grassland cover change better explained by single variable.The method actually controlled more variables (DEM, TEM, and PRE) through the division process.The results showed that the grassland cover change could be partly predicted by single variable, it was explainable by AD method.
In addition, the terrain factors used in this study (DEM, SLP, and ASP) were easy to measure with simple tools in practice, which accommodated the requirement of future study and grassland protection activities.
The analysis of the spatial-temporal changes of the NTV in the TRHR derived new spatial and temporal characteristics that could support future studies and policy making.Some possible driving factors are the erosion in area AAB, human activity in areas BAB(I) and BBB(III), precipitation and erosion in area BAB(III), permafrost thaw, rodent and "black soil beach" in area CAB, warm-drying tendencies in area CBB, and others, however further studies are required to strengthen the correlation of the relationships.The NTVV showed inverse tendencies towards the same factor in some areas.For example, the trend in areas BBB(III) and CAB to slope.The NTVV also showed opposite tendencies on both sides of the peak values in some areas.For example, there was an extreme value at 4500 m and the NTVV illustrated inverse trends on both sides of the value in area BAB(I).These phenomena require further study, and in addition, the NTVV rules and driving factors from more study areas, which were not deeply discussed in this study, need further investigation.More natural factors could be considered during the study AD process, such as evapotranspiration, land cover, and other related factors.
In a sustainability based ecological security perspective, human activities such as agriculture, livestock grazing, and infrastructure construction have resulted in the loss of natural ecosystems and the degradation of ecosystem services [57][58][59].The relationship between grassland cover change and terrain factors derived in this study might be helpful to indicate the degradation areas in TRHR.In construction, farming, and stocking activity planning in future, the fragile areas could be excluded, while natural reserves could be established in these areas.
Meanwhile, the second stage of the Conservation and Restoration Project has been published and applied concurrently [20], with clearly planned goals for the protection of glaciers, the management of "black soil beach", and the suppression of pests and rodents.It also guaranteed rotation grazing and the return of grazing land to grassland, moreover, it extended the eastern border of the protection zone.It could be predicted that in the future the overall NTV in the TRHR will increase, nevertheless, there will be constant ecological stresses due to high elevation, permafrost thaw, grassland erosion, and other environmental factors [60,61].

Conclusions
This study derived the NTV change characteristics in the TRHR via the AD and TA algorithms for a 17-year study period.Then the CART, RF, BAYE, and SVM machine learning methods were applied to evaluate the AD and TA results.The TRHR was divided into 12 study areas based on DEM, TEM, and PRE.The overall NTV rose 0.91% during the 17-year study period.In five-of-12 divided areas, the NTV demonstrated decreasing tendencies, while increases were observed in other areas.We focused on the severely decreased and variable change areas, and the results illustrated that the NTV showed apparent trends with regards to changes in ASP, DEM, or SLP.
Through ML evaluation, the AD and TA method applied in this study significantly increased the interpretation of grassland cover changes over the independent variables of DEM, SLP, and ASP.All the variables were convenient to measure through simple tools in practice.
In comparison with former study with AD methods based on administrative areas or basins in large or high-altitude difference regions, the division method based on natural conditions lead to better results for deriving the grassland cover change trends by single terrain factor together with TA method.In future research of grassland cover on large or high elevation difference areas, the methods applied in this study may help to provide better divisions and obtain improved results.
Simultaneously, the grassland cover changes of the TRHR in the QTP are complicated but unstable due to large areas, high elevations, low precipitations, thin soil layer, and other factors.Therefore, the influence of subsequent warm-wetting processes on vegetation will be difficult to predict, and thus, the continuous monitoring of grassland cover will be necessary.

Figure 1 .
Figure 1.The map of the location of the Three-River Headwaters Region (TRHR), the distribution of elevation in TRHR, and the map of Percent Non Tree Vegetation variation(NTVV) of MOD44B in TRHR.

Figure 1 .
Figure 1.The map of the location of the Three-River Headwaters Region (TRHR), the distribution of elevation in TRHR, and the map of Percent Non Tree Vegetation variation(NTVV) of MOD44B in TRHR.

Figure 2 .
Figure 2. Map of the low elevation, low temperature, low precipitation area (AAA); the low elevation, low temperature, high precipitation area (AAB); the low elevation, high temperature, low precipitation area (ABA); the low elevation, high temperature, high precipitation area (ABB); the moderate elevation, low temperature, low precipitation area (BAA); the moderate elevation, low temperature, high precipitation area (BAB); the moderate elevation, high temperature, low precipitation area (BBA); the moderate elevation, high temperature, high precipitation area (BBB); the high elevation, low temperature, low precipitation area (CAA); the high elevation, low temperature, high precipitation area (CAB); the high elevation, high temperature, low precipitation area (CBA); and the high elevation, high temperature, and high precipitation area (CBB).

Figure 3 .
Figure 3. (a) the counted raster data; (b) the data that determines the counting region; and (c) the result of the mean value of the TA algorithm.

Figure 3 .
Figure 3. (a) the counted raster data; (b) the data that determines the counting region; and (c) the result of the mean value of the TA algorithm.

Figure 4 .
Figure 4. Specific study procedures of the AD, TA, and ML analysis methods.

Figure 4 .
Figure 4. Specific study procedures of the AD, TA, and ML analysis methods.

Figure 6 .
Figure 6.(a) Map of the NTVV in the AAB area during the 17-year study period; and (b) the correlation between the NTVV and aspect in area AAB.

Figure 6 .
Figure 6.(a) Map of the NTVV in the AAB area during the 17-year study period; and (b) the correlation between the NTVV and aspect in area AAB.

Figure 7 .
Figure 7. (a) Map of the NTVV in area BAB(I) during the 17-year study period; (b) the correlation between the NTVV and the elevation in area BAB(I); (c) map of the NTVV in area BAB(II) during the 17-year study period; (d) the correlation between the NTVV and the elevation in area BAB(II); (e) map of the NTVV in area BAB(III) during the 17-year study period; and (f) the correlation between the NTVV and slope in area BAB(III).

Figure 7 .
Figure 7. (a) Map of the NTVV in area BAB(I) during the 17-year study period; (b) the correlation between the NTVV and the elevation in area BAB(I); (c) map of the NTVV in area BAB(II) during the 17-year study period; (d) the correlation between the NTVV and the elevation in area BAB(II); (e) map of the NTVV in area BAB(III) during the 17-year study period; and (f) the correlation between the NTVV and slope in area BAB(III).

Figure 9 .
Figure 9. (a) Map of the NTVV in area CAB during the 17-year study period; and (b) the correlation between the NTVV and slope in area CAB.

Figure 8 .
Figure 8.(a) Map of the NTVV in area BBB(I) during the 17-year study period; (b) the correlation between the NTVV and the elevation in area BBB(I); (c) map of the NTVV in area BBB(II) during the 17-year study period; (d) the correlation between the NTVV and the elevation in area BBB(II); (e) map of the NTVV in area BBB(III) during the 17-year study period; and (f) the correlation between the NTVV and slope in area BBB(III).

Figure 8 .
Figure 8.(a) Map of the NTVV in area BBB(I) during the 17-year study period; (b) the correlation between the NTVV and the elevation in area BBB(I); (c) map of the NTVV in area BBB(II) during the 17-year study period; (d) the correlation between the NTVV and the elevation in area BBB(II); (e) map of the NTVV in area BBB(III) during the 17-year study period; and (f) the correlation between the NTVV and slope in area BBB(III).

Figure 9 .
Figure 9. (a) Map of the NTVV in area CAB during the 17-year study period; and (b) the correlation between the NTVV and slope in area CAB.

Figure 9 .
Figure 9. (a) Map of the NTVV in area CAB during the 17-year study period; and (b) the correlation between the NTVV and slope in area CAB.

Figure 10 .
Figure 10.(a) Map of the NTVV in area CBB during the 17-year study period; and (b) the correlation between the NTVV and slope in area CBB.

Figure 10 .
Figure 10.(a) Map of the NTVV in area CBB during the 17-year study period; and (b) the correlation between the NTVV and slope in area CBB.

Table 1 .
Categorized research AD in the TRHR.
investigated the NTV variations (NTVV) and correlate drivers in areas with similar natural conditions.

Table 1 .
Categorized research AD in the TRHR.

Table 2 .
Research AD categories in the TRHR.

Study Areas * Code Combination Definitions 1st Letter: DEM (m) 2nd Letter: TEM ( • C) 3rd Letter: PRE (mm)
The areas and the NTVV over the 17-year study period are presented in Table3. *

Table 3 .
Research AD details in the TRHR.

Table 4 .
Details of ML analysis in further studied areas and TRHR.

Table 5 .
Details of RF ML results in further studied areas with (W) and without (WO) TA method.