Land use and land cover change detection and spatial distribution on the Tibetan Plateau

Regarded as the third pole of the Earth, the Tibetan Plateau (TP) is a region with complex terrain. Vegetation is widely distributed in the southeastern part of the plateau. However, the land use and land cover changes (LULCC) on the TP have not been sufficiently studied. In this study, we propose a method of studying the dynamic changes in the land cover on the TP. Landsat OLI images (2013 and 2015) were selected to extract the LULCC information of Nyingchi County, the DEM was used to extract objects’ land curved surface area and analyze their three-dimensional dynamic change information, which realized a four-dimensional monitoring of the forestry information on time and spatial level. The results showed that the forest area in 2015 decreased by 7.25%, of which the coniferous forest areas decreased by 25.14%, broad-leaved forest areas increased by 12.65%, and shrubbery areas increased by 14.62%. Compared with traditional LULCC detection methods, the change detection is no longer focused on the two-dimensional space, which helps determine the three-dimensional land use and land cover changes and their distribution. Thus, dynamic spatial changes can be observed. This study provides scientific support for the vegetation restoration and natural resource management on the TP.

types are the foundation for the study of vegetation coverage and dynamic changes. Therefore, studying the alpine vegetation and its dynamic changes is important to understanding the alpine ecological environment and climate change.
Studying the relationship between the vegetation distribution and the terrain condition could provide guidance for the exploration of the heterogeneity of the forest land and spatial distribution of the forestry resources. A digital elevation model (DEM) is a quantitative representation of the Earth's surface and the spatial information data, so it has become a vital source of information about the elevation, slope, and terrain in scientific investigations 27,28 . The remote sensing technique and DEM data have been widely used in land feature investigations, and these techniques and data can help improve the analytical ability of geoscience research 29 . For the past few decades, the TP has suffered from climate change and increasing temperatures, which have been aggravated in recent years. However, only a few studies have been conducted on the LULCCs on the TP, and precise and adequate information about land use change detection is extremely important for preferable management 30 . Therefore, in this study, a basic method for land cover change detection was applied to extract the land cover change information from 2013 to 2015, and Landsat OLI images were used as the main satellite data source. Previous studies rarely consider the impact of the topography on the land use type changes monitoring, the digital elevation model (DEM) is rarely used to study and extract the objects' land curved surface area. Therefore, the three-dimensional spatial information of the forestry resources changes has not been finely deep mined. Thus, in this study, the DEM data were used to analyze the spatial distribution of and dynamic changes in the land cover on the TP. The objective of this study was to accurately describe the land cover changes in the alpine vegetation on south-eastern TP from 2013 to 2015. As the terrain factor is considered in this study, the land use type and land cover change information are not confined in the two-dimensional plane level, which realized a four-dimensional monitoring of the mountainous forestry information, on the time and spatial level. With combination use of DEM and Landsat OLI images, the three-dimensional terrain visualization is realized, a real three-dimensional forestry land cover distribution and change monitoring were generated. In this study, the mapping relationship between the forest vegetation and the image features was established, which can improve the accuracy of information's spatial expression, it also can help the forest resource monitoring in the region of complex terrain. The results of this study are significant for natural resource conservation and vegetation resource management in alpine areas.

Results
Nyingchi County (93°27′-95°17′ E, 29°21′-30°15′ N) is located in the southeastern part of the Tibetan Autonomous Region of China (Fig. 1). The Yarlung Zangbo River, which is the largest river in China and one of the highest rivers in the world, flows through Nyingchi County. The study area is in a vegetation covered high mountain area with large topographic relief. The average altitude of the study area is 3000 m. The highest summit of the Nyingchi County is Gyala Peri Mountain, which is about 7294 m above sea level. It is the 85 th highest mountain in the world. The southwestern warm and humid air makes Nyingchi County a fertile area with abundant species.
Land use and land cover change detection. Based on the forestry inventory data and the field sampling data, the land cover types in the study area were divided into five types: coniferous forest (CF), broadleaf forest (BF), shrubbery (SH), non-forestry land (NF), and water (WA). Using the ENVI 5.3 platform, the  Maximum likelihood classification method was applied to extract the land cover information of the Landsat OLI images in 2013 and 2015, the forestry resource inventory data and sampling data were used as classification and validation samples, the confusion matrix was used to verify the classification accuracies, according to the verified results, the classification accuracies were 85.79% (2013) and 87.17% (2015). After obtaining the land cover classification results for the study area, the change detection statistics were calculated to obtain the land use transition matrix of the study area. Which is shown in Table 3. In Table 1, the diagonal elements are the plane projection areas of the land use types that did not change, i.e., land use unchanged (LUU), the off-diagonal elements are the plane projection areas of the land use types that did change, i.e., land use changed (LUC).
Dynamic land cover changes on the curved surfaces. Based on the curved land surface areas of the 20 LUCs and 5 LUUs, a new land use transition matrix was obtained. The statistical land use change data for 2013 and 2015 is shown in Table 2, and the statistical land use area data for 2013 and 2015 is shown in Table 3.
According to the classification results, over 80% of the land in the study area is covered by forests. Based on the classification results for 2015, 37.46% of the study area's land was covered by CF, which is the most widely distributed land use cover type. The SH accounted for 26.20% of the land area, and 22.36% of the land in the study area was BF. The rest of the land area was NF, accounting for about 12.53%, and only 1.54% of the land in   The statistical data show that the total forestry area decreased by 7.25%, of which the SH area increased by 14.62% and the CF area decreased by 25.14%. According to the comparison, the NF had a significant increment. Its area increased by 104.71% in 2015, while the WA increased by 31.29%. Table 3 shows the land cover area transfer data for the study area for 2013 and 2015. Based on the data presented in Table 3 Within the 4177-4576 m altitude gradient, the CF area's reduction level decreased as the altitude increased. Most of the increases in the CF area occurred in the 4577-4876 m altitude gradients. In the 3177-3976 m altitude gradients, the BF area increased, but the increase was only slightly correlated with the altitude. However, as the altitude increased to 3977-4376 m, the BF area decreased. The SH area increased in the 3777-4476 m altitude gradient, whereas in the 4477-4876 m altitude gradient, it decreased. The SH was mainly distributed at lower altitudes; however, very little SH occurred in the 3277-3776 m altitude gradient. Most of the increases in the NF area occurred in the 2977-3276 m and 4477-4876 m altitude gradient; very little NF occurred in the 3277-4376 m altitude gradient. The WA area increased in the 2977-3176 m altitude gradient, due to the interception of the dam and water storage project, there is a significant increase in the area of upstream water bodies. Most of the water body increase areas were transferred from the CA and SH, not too much BF can be found transferred into WA in this altitude gradient; very little WA occurred in the 3177-5061 m altitude gradient. Overall, the areas change of WA is closely related to human activity at the lower altitude region, which include socio-economic construction and infrastructure construction.
Based on LUC and LUU DEM data for the CF, BF, and SH, the sum of the relevant curved land surface areas in the 21 altitude gradients were calculated. Based on the line charts, their distributions were analyzed, and the cloud point data were used to create a three-dimensional topographic map. The three-dimensional topographic map of the CF, BF, and SH in both the LUC and LUU areas overlap with the three-dimensional topographic map of the study area. Therefore, the dynamic changes in the forest land on three-dimensional map were obtained. Based on this three-dimensional map, the spatial dynamic changes in the forest land cover were analyzed.
According to Fig. 4, in the 2977-4176 m altitude gradient, as the altitude level increased, the unchanged CF area decreased and very few unchanged CF areas occurred in the 4777-5061 m altitude gradient. The water vapor carried by the monsoon crosses the mountain and heats up on the leeward slope, there is more precipitation in the high slopes than in the valleys. Therefore, the CF conservation areas are widely distributed at high altitudes. In the 3077-3976 m altitude gradient, the CF transformed into BF occurred in the 3077-3976 m altitude gradient. The CF mainly transformed into SH in the 3077-3176 m and 3877-4576 m altitude gradient; the CF mainly transformed into NF in the 2977-3176 m altitude gradient; and the CF only transformed into WA in the 2977-3176 m altitude gradient. Figure 5a shows that the unchanged CF area was mainly distributed on the northern aspect of the mountains, and very little of the unchanged CF area was distributed in the southern and eastern aspect of the mountains. Figure 5b shows that the CF mainly transformed into BF in the western part of the study area, and its distribution has no direct relationship with mountains. Thus, the CF transformed into BF for different directions in the mountains, and it was mainly distributed in the 3200-3800 m altitude gradient. In the northern part of the study area, the CF was transformed into SH, mainly on the eastern and southern aspects of the mountain.   www.nature.com/scientificreports/ Along the river, the CF transformed into NF. Very little CF transformed into WA, and this transformation was mainly caused by water storage projects. Figure 6 shows that the transformation of BF into CF was widely distributed in the 3077-4276 m altitude gradient. Above 4277 m, this change rarely occurred. The unchanged BF area was mainly distributed in the 3077-4176 m altitude gradient. The transformation of BF into SH was widely distributed in the 3577-4476 m   Figure 7a shows that the unchanged BF area was mainly distributed in the western part of the study area, and most of the unchanged BF area was located on the eastern, western and northern aspects of the mountains. Very little of the unchanged BF area occurred on the southern aspect of the mountains. Figure 7b shows that the transformation of BF into CF mainly occurred in the central and eastern parts of the study area, at altitudes of 3000 to 4000 m. The distribution was not directly related to the aspect direction. The transformation of BF into SH mainly occurred in the northwestern part of the study area, and very little occurred in the southern and eastern parts of the study area. Barely any of this transformation occurred on the western and northern aspects of the mountains. The transformation of BF into NF mainly occurred along the riverbank. As with the CF, very little BF was transformed into WA in the study area. Figure 8 shows that above 4077 m, the SH was transformed into both CF and BF, and there were also SH areas that remain unchanged. In the 2977-3276 m and 4277-4976 m altitude gradients, the SH was transformed into NF, and the SH was only transformed into WA in the 2977-3176 m altitude gradient. Figure 9a shows that the unchanged SH area was mainly distributed in the northern and southwestern parts of the study area, and they mainly occurred on the western and southern aspects of the mountains. Figure 9b shows that the SH was mainly transformed into CF on the northern and western aspects of the mountains. The SH was widely transformed into CF in the study area, but the transformation was scattered across the study area, and no large areas of transformation occurred. Both at high and low altitudes, the SH was transformed into NF, and mainly occurred on the southern and eastern aspects of the mountains. The SH was transformed into WA where water storage projects were implemented.

Discussion
In this study, Nyingchi County was selected as the study area. The Landsat OLI images and DEM datasets were used as the main datasets for studying the LULCCs in Nyingchi County. The main goal of this study was to determine the spatial distributions of and changes in the land cover types in the study area. The results provide a better understanding of the dynamic LULCCs that have occurred in recent years. In this study, the land cover change detection was no longer restricted to the planar scale. As the terrain conditions of the TP are quite unique compared to other places on Earth, in order to study the spatial changes in the land cover types in the study area, the use of DEM data was introduced in this study. Twenty LUC areas and five LUU land use areas on curved land surface areas were calculated. The forest land cover change results were optimized, and the dynamic changes in the forest land cover areas were analyzed. Based on the result, the total forest area in the study area in 2015 decreased. In 2015, the CF was mainly transformed into SH and BF, while the BF was the main land cover that transformed into CF. The total CF area decreased by 25.14% compared with the CF areas in 2013. For the SH in the study area, most of the SH was transformed into NF in 2015, and most of the areas of SH increase in 2015 were transformed from CF. The SH area increased by 14.62% in 2015. In 2015, the BF area increased by 12.65%, of which, most of the BF area was transformed into CF, and most of the BF increased were due to the transformation of CF. The area of CF transformed into BF was about twice that of the area of BF transformed into CF.
The highest elevation in the study area is 5061 m. Based on altitude, the study area was divided into 21 altitude gradients. The curved land surface areas of the land cover types in the different altitude gradients and the changed and unchanged areas were calculated. Based on the results of the change detection, the spatial dynamic changes in the land cover types were analyzed. The results indicate that the transformation of BF into CF mainly  Most of the previously studies on the monitoring the dynamic changes of the forestry land cover type focused on two-dimensional level, which cannot accurately reveal the spatial distribution of vegetation coverage. Due to terrain impacts on the vegetation distribution, for the forest that located at the mountainous region, for study of land use classification and change detection, the spectral information of the remote sensing images was not sufficient. Therefore, the DEM data was used in the study to extraction the spatial distribution and dynamic change of the different land use type of the study area in 2013 and 2015. The results of this study reveal that the proposed method is an effective means of studying the LULC in complicated terrains. These results provide a better understanding of the spatial distributions of and changes in the land cover types, which provides a more intuitive understanding of the dynamic changes in land cover type on the TP. According to the study results, monitoring the mountainous forest on spatio-temporal level and deep mining the three-dimensional information of the forest resource dynamic changes can help formulate accurate forest protection policies and realize the sustainable management of the forest resources. However, there are still several problems regarding the dynamic change in the land cover type that should be studied further. In this study, we focused on a short time series of LULC, and the driving factors of the land use change were not investigated and discussed in this study. Further research on the driving factors and the LULC's impact on the environment and climate is necessary.

Conclusions
In this study, the LULCCs from 2013 to 2015 in Nyingchi County on China's Tibetan Plateau were analyzed. In this study, DEM data were used in the calculations of the changed and unchanged curved land surface areas, and thus, the land cover change results were modified. In the places like the TP, where the terrain is complex, the calculation of the curved land surface areas is close to the real conditions. The altitude of the study area was divided into 21 gradients, and the land cover changes from 2013 to 2015 were analyzed. The method proposed in this study can effectively extract the spatial dynamic changes in the land cover types in the study area. Based on the three-dimensional graphics of the land cover change results, the graphics reveal the land cover change information for the true three-dimensional terrain conditions. The three-dimensional graphics obtained using the method proposed and applied in this study have a rather good visual effect. Because there are no shaded sides, the distributions of both the changed and unchanged areas can be easily observed, and thus, these graphics can be effectively used to analyze the spatial dynamic land cover changes.

Materials and methods
Satellite images and data processing. Landsat 8 images were used in this study. The Landsat 8 satellite was launched by the United States Geological Survey (USGS) and the National Aeronautics and Space Administration in 2013, and it contains two sensors: the operational land imager (OLI) and the thermal infrared sensor. The Landsat OLI image contains nine bands with a resolution of 30 m, and a full-color band with a resolution of 15 m. The images used in this study were downloaded from the USGS and the Internal Scientific Data Service Platform of China. Based on the study's goal of land cover change detection, two sets of images taken in 2013 and 2015, with cloud cover of less than 5%. The Landsat OLI images were obtained during the vegetation growing season. Before the remote sensing image classification and land cover change detection, the images were processed using ENVI 5.3 software. The image processing included radiometric, atmospheric, geometric, and topographic corrections. After the atmospheric correction was made, the vegetation spectrum curve was close to the real vegetation spectrum. Images that cover the Yarlung Zangbo River Basin were obtained.
The DEM dataset was also used in this study. The spatial resolution of the DEM was 1 arcsec. In order to make it consistent with the remote sensing images, the DEM dataset was resampled. The resample spatial resolution   www.nature.com/scientificreports/ was 30 m, and cubic convolution interpolation was used for the dataset resampling. However, even though the remote sensing images were geo-referenced, the images and the DEM dataset did not line up exactly. In order to match the different images, the remote sensing images were resampled based on the DEM dataset using the  where y i (i = 1, 2, …, n) represents the abscissa values and x i (i = 1, 2, …, n) represents the ordinate values of the spatial rectangular coordinate system O-XYZ. h k (k = 1, 2, …, n*m) is the height of the altimetric point, n is the altimetric points position on the east--west line; and m is the altimetric points position on the north-south line. The corresponding coordinate values can be calculated using m and n, which are obtained based on the DEM matrix.
The spatial resolution of the DEM is 30 m, which means y i is 30 m larger than y i−1 , and x j is 30 m larger than x j−1 . Thus, the complete data grid is defined by the following matrixes: (2) y 1×n = y 1 y 2 · · · y n (3) x 1×n = x 1 x 2 · · · x m (4) y m×n =      y 1 y 2 · · · y n y 1 y 2 · · · y n . . . . . . . . . . . . www.nature.com/scientificreports/ Each grid point is related to a height value, and the height value matrix is defined as The data set y, x, h m×n can be expressed as the function h = f y, x . In this study, the gradients in both the east-west direction and the south--north direction were calculated using the gradient function fy, fx = gradient h, step, step . The average space in every direction is a step, and because the resolution of the DEM data is 30 m, the step value was set as 30 before the calculation. The element of the curved surface areas can be calculated by formula (Eq. 7). The sum of all the elements is the area of the curved land surface. The curved land surface areas of 20 LUCs and 5 LUUs were used to analyze the dynamic change in the land cover in the study area.