Vegetation-Ice-Bare Land Cover Conversion in the Oceanic Glacial Region of Tibet Based on Multiple Machine Learning Classiﬁcations

: Oceanic glaciers are one of the most sensitive indicators of climate change. However, remotely sensed evidence of land cover change in the oceanic glacial region is still limited due to the cloudy weather during the growing season. In addition, the performance of common machine learning classiﬁcation algorithms is also worth testing in this cloudy, frigid and mountainous region. In this study, three algorithms, namely, the random forest, back-propagation neural network (BPNN) and convolutional neural network algorithms, were compared in their interpretation of the land cover change in south-eastern Tibet and resulted in three ﬁndings. (1) The BPNN achieves the highest overall accuracy and Kappa coe ﬃ cient compared with the other two algorithms. The overall accuracy was 97.82%, 98.07%, 98.92%, and 94.63% in 1990, 2000, 2007, and 2016, and the Kappa coe ﬃ cient was 0.958, 0.959, 0.980, and 0.918 in these four years, respectively. (2) From 1990 to 2000, the dominant land cover was ice at the landscape level. The landscape fragmentation decreased and the landscape aggregation increased. From 2000 to 2016, the dominant land cover transformed from ice to vegetation. The vegetation aggregation increased, while the ice aggregation decreased. (3) When the elevation was less than 4 km, the vegetation was usually transformed into bare land; otherwise, the probability of direct transformation between vegetation and ice increased. The ﬁndings on the land cover transformation in the oceanic glacial region by multiple classiﬁcation algorithms can provide both long-term evidence and methodological indications to understand the recent environmental change in the “third pole”.


Introduction
As an important indicator of global warming, the long-term change in glaciers can provide strong evidence for the understanding of the recent environmental changes [1,2]. With the warming climate and the increasing intensity of human activities, most glaciers around the world have been receding [3]. Since 1960, global glaciers have experienced a strong contraction in the past 50 years, with a contraction rate of 11.3%; the annual percentage change was 0.35% [4]. The fast and continuous shrinkage of glaciers has some benefits on the environment and human activities, including vegetation growth, economic development and increasing the water and energy supply [5,6]. However, many hazards have also resulted from glaciers, such as ice avalanches, debris flows and glacial lake outburst floods, which

Study Area
The study area is a part of Bomê County, Tibet, China ( Figure 1). Bomê County is located in south-eastern Tibet. Its geographical coordinates are between 29 • N and 30 • N and between 94 • E and 96 • E. It is characterized by a subtropical mountain humid monsoon climate, with an annual mean temperature of 8.5 • C and an annual precipitation of 900 mm. The total land area of Bomê County is 16,795 km 2 , and the elevation ranges from 2139 m to 6462 m, with a height difference of more than 4000 m. Bomê County is one of the typical representative areas of oceanic glaciers in China and the oceanic glaciers in the area are extremely well developed. Due to its special geographical location and climatic characteristics, glaciers easily accumulate and melt, which increases the difficulty of accurate temporal observations.

Data and Samples
The study area is fully covered by Landsat path 135 and row 39. Four Landsat image scenes of the different periods (Table 1) were loaded from the Geospatial Data Cloud (http://www.gscloud.cn/) and are all "Tier 1" products, representing the highest quality imagery meeting certain radiometric and geometric requirements.  The four images of the study area, with size 5777 × 4013 pixels, are in the UTM coordinate system geo-referenced to WGS84. We resampled all the image data to a 30 m resolution to ensure that all bands matched the spatial resolution of both the training and validation data. We input all the bands into the classification network for calculation. In addition, the normalized difference vegetation index (NDVI), normalized difference snow index (NDSI), and normalized difference built-up index (NDBI) as input data are also put into the networks to increase the classification accuracy. The NDVI is one of most well-known and frequently used indexes computed from multi-spectral remote sensing [28]. The NDSI is defined as the difference between the reflectances of the green and short-wave infrared bands divided by the sum of the two reflectances [29]. The snow/ice has high reflectance in the visible light region and strong absorption in the SWIR region, which means that the NDSI can effectively distinguish snow from similar bright soils, vegetation and rocks, and clouds in TM images [30]. The NDBI was first proposed in 2003 to extract urban built-up land using Landsat images [31]. The index output ranges from -1 to +1, with positive values representing built-up areas [32]. The NDBI is unable The four images of the study area, with size 5777 × 4013 pixels, are in the UTM coordinate system geo-referenced to WGS84. We resampled all the image data to a 30 m resolution to ensure that all bands matched the spatial resolution of both the training and validation data. We input all the bands into the classification network for calculation. In addition, the normalized difference vegetation index (NDVI), normalized difference snow index (NDSI), and normalized difference built-up index (NDBI) as input data are also put into the networks to increase the classification accuracy. The NDVI is one of most well-known and frequently used indexes computed from multi-spectral remote sensing [28]. The NDSI is defined as the difference between the reflectances of the green and short-wave infrared bands divided by the sum of the two reflectances [29]. The snow/ice has high reflectance in the visible light region and strong absorption in the SWIR region, which means that the NDSI can effectively distinguish snow from similar bright soils, vegetation and rocks, and clouds in TM images [30]. The NDBI was first proposed in 2003 to extract urban built-up land using Landsat images [31]. The index output ranges from -1 to +1, with positive values representing built-up areas [32]. The NDBI is unable to separate urban areas from bare land, but we can still use it to extract bare land from other land cover types because built-up areas are rare in this region [33]. The NDVI (Equation (1)), NDSI (Equation (2)) and NDBI (Equation (3)) of each image were calculated using the Raster Calculator tool in ArcGIS 10.5: where B N is the reflectance of the near-infrared band (band 4 for the Landsat 5 TM and Landsat 7 ETM+ sensor and band 5 for the Landsat 8 OLI sensor); B R is the reflectance of the red band (band 3 for Landsat 5 TM/Landsat 7 ETM+ and band 4 for Landsat 8 OLI); B G is the reflectance of the green band (band 2 for Landsat 5 TM/Landsat 7 ETM+ and band 3 Landsat 8 OLI); B S is the reflectance of the short-wave infrared band (band 5 for the Landsat 5 TM and Landsat 7 ETM+ sensor and band 6 for the Landsat 8 OLI sensor). We used ArcGIS tools to randomly generate 20,000 random points in the study area (the 20,000 sample points on the four images were in the same position). We set the land cover types in this area into vegetation, ice, bare land, clouds, water, and shadows by setting thresholds for NDVI, NDSI, NDBI, and specific spectral bands. The judgment rules are as follows: (1) when a random point is NDVI> = 0.1 and NDVI> NDSI and NDVI> NDBI, it is judged as vegetation; (2) when a random point is NDSI> = 0 and NDSI> NDVI, It is judged as ice or cloud; (3) When a random point has NDBI> = 0.1 and NDBI> NDVI, it is judged as bare ground; (4) When a random point has NDVI <0 and blue band spectral value less than 120 (the blue band spectral value of the OLI sensor is greater than 10,000), it is judged as water; (5) Ice and clouds are distinguished by short-wave infrared. When it is less than 11,000 (for 2016 images only), it is judged as ice; otherwise, it is cloud. Since the other three images selected (1990, 2000 and 2007) were almost cloudless, visual interpretation was used to distinguish between clouds and snow. After the preliminary classification is completed, the four images (1990, 2000, 2007 and 2016) have 280, 567, 341, 1654 points left, respectively, which are inspected by manual judgment based on Google Earth high-resolution images. The number of training/validation points for different land cover types in four years is shown in Table 2, in which three-quarters of each category (75% of the total number of individual categories) were randomly selected as the training set, the remaining (25%) was used as the validation set. The training set and validation set of the CNN are square images of size 7 × 7 pixels with sample points as centre points.  T  V  T  V  T  V  T  V   Vegetation  4097  1365  2894  965  4216  1405  5657  1885  Ice  9218  3073  10198  3399  8588  2862  6090  2030  Bare land  1188  396  1431  477  1868  622  235  78  Cloud  37  12  63  21  10  3  2235  745  Water  53  17  86  28  33  11  84  28  Shadow  161  54  82  27  40  13  453  151 Remote Sens. 2020, 12, 999 6 of 20

Methods
The first machine learning algorithm we applied is random forest (RF), which contains multiple decision trees, and every decision tree is trained on different input data sets generated through a bagging procedure, ultimately creating an ensemble of independent experts. The final classification result is based on a collection of all of the trees [34]. To balance the complexity and computational efficiency of the network, the RF algorithm is set to contain 50 decision trees, and the depth of the trees is 5.
The second machine learning algorithm we applied is the BPNN, which consists of forward path and backward path. The forward path involves creating a feed-forward network, initializing the weights, simulation and training the network [17]. The feed-forward network has three parts, namely, input layers, hidden layers and output layers, and at the same time, the layers are directly connected by means of full connections. The backward path is mainly used for updating the network weights and biases that are learned to approximate the complex relationship between the input features and the output characteristics [17]. A specific cost function is used to minimize the difference between the real values and the predicted output [35]. In this study, the BPNN consisted of one hidden layer that had 12 computational nodes. The learning rate and the number of epochs were chosen as 0.01 and 80, respectively. We used a simple and commonly used mean square error loss function.
The third algorithm we used is the CNN algorithm, which was used to process and analyse large-scale sensor data or images by learning stationary characteristics at the local and global scales [23], with alternatively connecting convolutional layers and pooling layers to obtain the deep and abstract features. To balance the complexity and robustness of the network, the CNN was designed to consist of two repetitive parts and three full connected layers. There are two identical convolutional layers and a max pooling layer in each repetitive part, and each convolutional layer is followed by a rectified linear unit activation function. The input size of the CNN is 7 × 7 pixels, the filter size of the convolutional layer was set to 3 × 3 with one stride, and the number of convolutional filters was set to 16, 32, 64, and 128 in turn. The learning rate and batch size were chosen as 0.0001 and 800, respectively. In addition, the number of epochs was set to 200 to learn the features fully with back propagation.
To make full use of spectral information, the input data are the spectral values of all bands and the three indexes (NDVI, NDSI and NDBI). Thus, the input feature numbers of the four years are 10, 12, 10 and 14, respectively. The predicted output of the three networks is expected to be in six classes. Once the final model was trained well; it was applied to an independent set of validation data that was not used to train the models.

Accuracy Evaluation
We use producer's accuracy (PA), user's accuracy (UA), overall accuracy (OA) and Kappa coefficient to evaluate the performance of the three classifiers. The producer's accuracy and user's accuracy can be acquired from the confusion matrix. The producer's accuracy represents the proportion of samples that are correctly classified as category i among all the category i samples (a column of the confusion matrix). The user's accuracy indicates that among all the samples classified as category i, the measured type is indeed the proportion of category i samples. The overall classification accuracy indicates the proportion of samples that were correctly classified among all the samples. The overall accuracy is proportional to the area that is correctly mapped and is suitable for area estimation. The Kappa coefficient is a measure of the overall agreement in the matrix. In contrast to the overall accuracy, the Kappa coefficient also takes the non-diagonal elements into account and has become a widely used measure of classification accuracy [36]. Since the samples in this study were randomly selected from the entire study area, the classification accuracy can be better evaluated by selecting the above four accuracy evaluation indicators.

Landscape Indexes
Landscape indexes can effectively summarize the landscape pattern information within a certain range, reflect the structural composition and spatial configuration of different landscape types, and thus quantitatively analyse landscape characteristics and changing trends. The landscape indexes at the class level and landscape level for the four years (1990, 2000, 2007 and 2016) were calculated in Fragstats 4.2 to portray the characteristics of the spatial and temporal changes in the various landscape categories. The class-level indexes reflect the amount and spatial distribution of a single patch type, while the landscape-level indexes reflect the spatial pattern of the entire landscape mosaic [37]. In this study, the class-level indexes, including patch density (PD), mean patch area (AREA_MN), edge density (ED), landscape shape index (LSI), largest patch index (LPI) and aggregation index (AI), were used. In addition to the above indexes, the landscape-level indexes also include the contagion index (CONTAG) and Shannon's diversity index (SHDI). These indexes are widely used to represent landscape compositions and configurations, as shown in Table 3.
To a certain extent, AREA_MN and the PD both reveal the fragmentation of patches in the landscape elements. AREA_MN reflects the average patch size of the landscape elements. The smaller its value, the higher the fragmentation of the landscape. The PD indicates the fragmentation degree of landscape element patches. The larger the PD value is, the smaller the patch size per unit area and the higher the landscape heterogeneity and fragmentation degree.
The ED is mainly a description of the degree of segmentation of landscape unit composition. The greater the density of the landscape edges is, the greater the degree of landscape fragmentation by the edges. The LSI can describe the regularity and complexity of landscape patch boundaries and refers to the degree of difference between circles or squares with the same patch area and shape. The larger the index value is, the more irregular the shape of the patch, and the less likely it is to be disturbed by the outside world. The LPI quantifies the percentage of the total landscape area composed of the largest patches and is an indicator of the advantages of landscape types. Moreover, the impact of human activities on the landscape can be reflected to some extent by changes in this value. The AI reflects the degree of aggregation. When its value is 0, the degree of fragmentation is the largest. To some extent, it reflects the complexity of the landscape and reveals the degree of human interference in the landscape. CONTAG mainly describes the degree of convergence or extension trend in different patch types in the landscape. In general, a high spread value indicates that a dominant patch type in the landscape has formed good connectivity; otherwise, it indicates that the landscape is more fragmented. The SHDI is the quantifying diversity at the landscape level. When the landscape has only one patch, the SHDI equals 0, and when the patch type increases, the value of the SHDI also increases. The evolution process of landscape patterns can be determined by investigating the basic characteristics, morphological changes and spatiotemporal evolution of the various landscape categories [38]. Table 3. Landscape pattern indexes used in this study.
Aggregation Index (AI) g ii = number of like adjacencies (joins) between pixels of patch type (class) i based on the single-count method.
Max->g ii = maximum number of like adjacencies (joins) between pixels of patch type (class) i (see below) based on the single-count method. P i = proportion of landscape comprised of patch type (class) i.

Contagion Index (CONTAG) [39]
P i = proportion of the landscape occupied by patch type (class) i. g ik = number of adjacencies (joins) between pixels of patch types (classes) i and k based on the double-count method. m = number of patch types (classes) present in the landscape, including the landscape border if present.

Classification Accuracy
The classification results are shown in Figure 2. Due to the small number of samples in the categories of water and shadows, only the vegetation, ice, bare land, and clouds are analysed for accuracy here. A summary of the results for three tested models considering four classes is shown ( , the BPNN improves the OA by more than 1.36%/0.29%/1.06%/2.08% and 1.99%/1.12%/2.11%/1.42%, and the Kappa coefficients by more than 0.027/0.004/0.019/0.034 and 0.039/0.024/0.038/0.021 over the RF model and the CNN model for the four years, respectively. Although the BPNN algorithm is superior to the other two algorithms in terms of the OA and Kappa coefficients, it did not achieve the best accuracy in terms of the PA and UA. For example, for a single kind of vegetation, the UA values of the BPNN algorithm for 1990, 2000, 2007, and 2016 were 97.96%, 97.67%, 99.36%, and 96.60%, respectively, which were higher than those of the RF algorithm (97.73%, 95.07%, 99.04%, and 90.95%, respectively) and that of the CNN algorithm (97.13%, 96.55%, 98.95%, and 96.71%, respectively). For ice, the UA values of the BPNN algorithm on the images at four different periods were 98.65%, 99.44%, 99.48%, and 95.15%, which were 1.01%, 0.06%, 0.59%, and 0.59% higher than the UA values of the RF algorithm, respectively, and were 1.57%, 0.84%, 1.1%, and 1.21% higher than the UA values, respectively, of the CNN algorithm. However, for 2000 and 2016, the PAs of vegetation were 95.75% and 96.45%, which were slightly lower than the RFs of 95.95% and 97.61%, respectively. For ice, the PA value was only 98.67% in 2016, which was higher than the 97.68% PA from the RF model. Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 20

Changes in the Landscape Pattern
To further analyse the landscape pattern characteristics of the land cover types in the study area, we adopted interpolation to remove clouds and shadows. First, the BPNN algorithm classification map with the highest classification accuracy is selected as the base map, and the RF and CNN algorithm classification maps are used as the first and second auxiliary maps, respectively. Second, by overlaying the base map with the first auxiliary map, most of the cloud and shadow areas in the base map are replaced by the land cover types in the RF algorithm classification map, and most of the remaining areas are replaced by the land cover types in the second auxiliary map. Third, the maps from adjacent years are superimposed to remove the remaining cloud and shadow types by temporal replacement. The third step would be performed only if the base map does not completely remove the clouds and shadows after the second step. By interpolating the classification map, the influence of clouds and shadows on the analysis of landscape features is eliminated so that the landscape features of the land cover map can be effectively analysed (Figure 3). Sections 3.2.1 and 3.2.2 analyse the land cover characteristics of the oceanic glacial regions from the class level and landscape level, respectively.
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 20 To further analyse the landscape pattern characteristics of the land cover types in the study area, we adopted interpolation to remove clouds and shadows. First, the BPNN algorithm classification map with the highest classification accuracy is selected as the base map, and the RF and CNN algorithm classification maps are used as the first and second auxiliary maps, respectively. Second, by overlaying the base map with the first auxiliary map, most of the cloud and shadow areas in the base map are replaced by the land cover types in the RF algorithm classification map, and most of the remaining areas are replaced by the land cover types in the second auxiliary map. Third, the maps from adjacent years are superimposed to remove the remaining cloud and shadow types by temporal replacement. The third step would be performed only if the base map does not completely remove the clouds and shadows after the second step. By interpolating the classification map, the influence of clouds and shadows on the analysis of landscape features is eliminated so that the landscape features of the land cover map can be effectively analysed (Figure 3). Sections 3.2.1 and 3.2.2 analyse the land cover characteristics of the oceanic glacial regions from the class level and landscape level, respectively.   At the landscape class level, the LPI and AI are used here to reflect the degree of aggregation of the landscape elements. The LPI and AI values for ice both increased from 1990, reached the highest peak in 2000, and declined after 2000. In contrast, the LPI and AI values of vegetation both decreased from 1990. They reached the lowest points in 2000 and increased after 2000. Additionally, the AI value for ice was greater than 95 for all four years. Compared with ice, the AI value for vegetation for the four years was lower, but all were greater than 90. Moreover, with respect to the change in the AI values for ice and vegetation, there was an opposite trend. It can be inferred that the degree of ice and vegetation accumulation was high, but the degree of ice accumulation was much higher than that of vegetation. In 2000, the ice had the highest concentration, while the concentration of vegetation was the lowest. From 2000 to 2007, the LPI for ice was much higher than that for vegetation, but the LPI for vegetation exceeded that for ice in 2016. It can be guessed that from 1990 to 2007, ice was the At the landscape class level, the LPI and AI are used here to reflect the degree of aggregation of the landscape elements. The LPI and AI values for ice both increased from 1990, reached the highest peak in 2000, and declined after 2000. In contrast, the LPI and AI values of vegetation both decreased from 1990. They reached the lowest points in 2000 and increased after 2000. Additionally, the AI value for ice was greater than 95 for all four years. Compared with ice, the AI value for vegetation for the four years was lower, but all were greater than 90. Moreover, with respect to the change in the AI values for ice and vegetation, there was an opposite trend. It can be inferred that the degree of ice and vegetation accumulation was high, but the degree of ice accumulation was much higher than that of vegetation. In 2000, the ice had the highest concentration, while the concentration of vegetation was the lowest. From 2000 to 2007, the LPI for ice was much higher than that for vegetation, but the LPI for vegetation exceeded that for ice in 2016. It can be guessed that from 1990 to 2007, ice was the main dominant land cover category in the study area, and in 2016, the vegetation area exceeded that of ice.

The Landscape-Level Pattern
The AREA_MN and PD can reflect the degree of fragmentation in a landscape. As shown in Figure 5, on the one hand, AREA_MN was approximately 5.0 in 1990, but it rose sharply to 6.8 in 2000 and then dropped to approximately 5.0 in 2016. On the other hand, the PD in 1990 was greater than 6.5, dropped to 4.6 in 2000, and then returned to 6.5 in 2016. In 2000, AREA_MN and the PD reached the highest peak and the lowest valley, respectively, which indicates that landscape fragmentation was large in 2000. From 1990 to 2000, the landscape fragmentation in the study area decreased, and after 2000, the landscape fragmentation increased. The tendency of landscape fragmentation to decrease first and then increase coincides with the trend in the fragmentation of ice and vegetation. It was revealed that the fragmentation of ice and vegetation in the study area reached the minimum in 2000, and the fragmentation increased after 2000. Compared with 1990, the fragmentation did not change much in 2016.

The Transformation of Land Cover
As shown in Figure 6, the lower semicircle represents the three types of land cover in the previous year, the upper semicircle represents the three types of land cover in the next year, and the transfer between different colors represents the change in the type of land cover (see Figure 6a Veg-1990's "red" flow to Bare-2000 refers to the percentage of land cover types in 1990 as vegetation, and in 2000 the land cover type became bare land). The land cover structure of the study area greatly changed from 1990 to 2000 (Figure 6a). Among them, the change in ice was the most significant. Approximately 41% of the bare land and 10% of the vegetation were converted into ice, which greatly increased the ice coverage in the area. In addition, nearly 18.7% of the vegetation coverage area was converted into bare land. From 2000 to 2007, the change in ice was relatively small, and nearly half of the area of bare land was converted into vegetation coverage in 2007 (Figure 6b). Figure 6c shows the mutual transfer of land cover types in the study area of Bomê County from 2007 to 2016. We found that most of the bare land has been converted into vegetation, and the glaciers slightly decreased from 2007 to 2016. This finding implies that the glacier area continued to decrease, and the area of vegetation coverage increased significantly from 2000 to 2016. Figure 6d reveals that the land cover structure of Bomê County changed drastically from 1990 to 2016. For instance, nearly three-quarters of the bare land was converted into vegetation. Accordingly, 99.9% of the vegetation areas were unchanged, but the proportion of bare land self-transfer was very low. In summary, our results indicated that vegetation was restored after 2000, and the coverage of vegetation reached its highest The value range of the ED was between 30 and 36, which reflects that the segmentation of the patch edge was not high, and the segmentation change during 1990-2016 was not large. The LSI value varied between 80 and 94, showing that the landscape shape becomes more irregular and complex and that the landscape changed a little due to external disturbance. Therefore, from the changes in the LSI and ED values, we can reveal that compared with other years, the landscape type in this study area was more regular in shape, and the degree of fragmentation was the lowest in 2000. However, compared to 1990, the landscape type in 2016 was slightly simpler and less fragmented.
The LPI values show that the main landscape types increased from 1990 to 2000 and then decreased from 2000 to 2007, consistent with the changes in the LPI for the glaciers from 1990 to 2007, which also shows that ice is the dominant landscape in this period. In addition, the sharply increasing LPI value from 2007 to 2016 is likely to be affected by the type of vegetation because the LPI value for vegetation increased rapidly from 2007 to 2016, and vegetation became the dominant landscape type. From 1990 to 2016, the value of the landscape AI was not less than 94. In addition, according to Figure 4, the AI values for ice from 1990 to 2016 were all greater than 97, and the AI values for vegetation were all greater than 91, with a small fluctuation range. Therefore, according to the above analysis, it can be seen that the concentration of various landscape types in this study area is very high, and the dominant landscape types are ice and vegetation. Moreover, the dominance of ice was higher than that of vegetation from 1990 to 2007, and 2016 was just the opposite.
At the landscape level, the SHDI is calculated mainly because it represents the diversity in the entire landscape level. The SHDI values during 1990-2016 were relatively stable and exhibited a slight overall decrease without significant changes, indicating a slight decrease in landscape diversity. The CONTAG values for the four years in this study area exhibited certain fluctuations, with a maximum value of 70 and a minimum value of 64, indicating that the patch-like connectivity in the study area is good, but the differences vary from year to year. From 1990 to 2000, the value of CONTAG slowly increased to 68, and from 2000 to 2007, the CONTAG value decreased from 68 to 64 and quickly increased to 70 in 2016. The above shows that from 1990 to 2000 and from 2007 to 2016, the connectivity of patches in the study area continued to increase, and the connectivity was the smallest in 2000. In summary, since the dominant patch types in the study area are ice and vegetation, the connectivity for ice and vegetation varies from year to year, and the connectivity is the worst in 2007.

The Transformation of Land Cover
As shown in Figure 6, the lower semicircle represents the three types of land cover in the previous year, the upper semicircle represents the three types of land cover in the next year, and the transfer between different colors represents the change in the type of land cover (see Figure 6a Veg-1990's "red" flow to Bare-2000 refers to the percentage of land cover types in 1990 as vegetation, and in 2000 the land cover type became bare land). The land cover structure of the study area greatly changed from 1990 to 2000 (Figure 6a). Among them, the change in ice was the most significant. Approximately 41% of the bare land and 10% of the vegetation were converted into ice, which greatly increased the ice coverage in the area. In addition, nearly 18.7% of the vegetation coverage area was converted into bare land. From 2000 to 2007, the change in ice was relatively small, and nearly half of the area of bare land was converted into vegetation coverage in 2007 (Figure 6b). Figure 6c shows the mutual transfer of land cover types in the study area of Bomê County from 2007 to 2016. We found that most of the bare land has been converted into vegetation, and the glaciers slightly decreased from 2007 to 2016. This finding implies that the glacier area continued to decrease, and the area of vegetation coverage increased significantly from 2000 to 2016. Figure 6d reveals that the land cover structure of Bomê County changed drastically from 1990 to 2016. For instance, nearly three-quarters of the bare land was converted into vegetation. Accordingly, 99.9% of the vegetation areas were unchanged, but the proportion of bare land self-transfer was very low. In summary, our results indicated that vegetation was restored after 2000, and the coverage of vegetation reached its highest amount in 2016.

Relation between Land Cover Change and Vegetation
To study the mutual conversion relationship between vegetation, ice and bare land at different elevations, the conversion areas from 1990 to 2000, 2000 to 2007, 2007 to 2016, and 1990 to 2016 at different elevation ranges are plotted in Figure 7. Figure 7a shows that from 1990 to 2000, when the elevation was less than 4 km, the area converted from vegetation to ice was almost twice the area converted from bare land to ice. With increasing elevation, the area converted from bare land to ice increases, exceeding the area converted from vegetation to ice. From Figure 7b-d, when the elevation is less than 4 km, the area converted from bare land to vegetation is larger than the area converted from ice to vegetation. However, when the elevation is higher than 4 km, the area converted from ice to vegetation is larger than the area. We can speculate that as elevation increases, vegetation is more likely to convert directly into ice rather than bare land. At the same time, as the elevation increases, the type of land cover converted to vegetation is most likely to be ice rather than bare land. It can be concluded that with increasing elevation, the conversion probability between ice and vegetation is increasing.  Veg-1990 indicates that the type of land cover in 1990 was vegetation.

Relation between Land Cover Change and Vegetation
To study the mutual conversion relationship between vegetation, ice and bare land at different elevations, the conversion areas from 1990 to 2000, 2000 to 2007, 2007 to 2016, and 1990 to 2016 at different elevation ranges are plotted in Figure 7. Figure 7a shows that from 1990 to 2000, when the elevation was less than 4 km, the area converted from vegetation to ice was almost twice the area converted from bare land to ice. With increasing elevation, the area converted from bare land to ice increases, exceeding the area converted from vegetation to ice. From Figure 7b-d, when the elevation is less than 4 km, the area converted from bare land to vegetation is larger than the area converted from ice to vegetation. However, when the elevation is higher than 4 km, the area converted from ice to vegetation is larger than the area. We can speculate that as elevation increases, vegetation is more likely to convert directly into ice rather than bare land. At the same time, as the elevation increases, the type of land cover converted to vegetation is most likely to be ice rather than bare land. It can be concluded that with increasing elevation, the conversion probability between ice and vegetation is increasing.

Comparison of the Classification Algorithms
In this study, three machine learning algorithms were used to classify land cover types in Bomê County, a typical oceanic glacial region. The results show that the BPNN algorithm achieves better classification accuracy than the RF and CNN algorithms. The BPNN is a widely used feed-forward artificial neural network (ANN). By comparing the performances of RF and ANN for the Land use/cover classification, Ge et al. concluded that the highest OA was produced by the ANN (97.16%), which was higher than the RF (96.92%) [40]. Consistent with the classification results, the BPNN used in this study has a higher classification accuracy for the identification of glaciers in oceanic glacial regions. Therefore, the BPNN algorithm has great development potential in oceanic glacial regions. Yoo et al. made a detailed comparison between CNN and RF classifiers to find that CNN consistently resulted in higher accuracy than RF, and OA increased by 7%-18% [41]. However, the OA for the CNN algorithm is inferior to those of the BPNN and RF in this study. In the CNN-based scheme, the influence of neighbourhood pixels as input features is considered. However, a common practical problem based on deep learning classifiers is the class imbalance problem, which means that some neighbourhood classes have more examples in the training set than others. The overall accuracy measure is related to difficulties in the context of imbalanced data, and class imbalances can have an adverse effect on the training of deep learning classifiers [42]. This problem affects both the convergence of the training phase and the generalization ability of the model on the test set. In this study, due to the different neighbourhood sample sizes of the six categories, there is a certain class imbalance in the training and test samples in CNN compared with that in BPNN and RF, which may mislead the class label of training set and affect the CNN classification performance.

Comparison of the Classification Algorithms
In this study, three machine learning algorithms were used to classify land cover types in Bomê County, a typical oceanic glacial region. The results show that the BPNN algorithm achieves better classification accuracy than the RF and CNN algorithms. The BPNN is a widely used feed-forward artificial neural network (ANN). By comparing the performances of RF and ANN for the Land use/cover classification, Ge et al. concluded that the highest OA was produced by the ANN (97.16%), which was higher than the RF (96.92%) [40]. Consistent with the classification results, the BPNN used in this study has a higher classification accuracy for the identification of glaciers in oceanic glacial regions. Therefore, the BPNN algorithm has great development potential in oceanic glacial regions. Yoo et al. made a detailed comparison between CNN and RF classifiers to find that CNN consistently resulted in higher accuracy than RF, and OA increased by 7%-18% [41]. However, the OA for the CNN algorithm is inferior to those of the BPNN and RF in this study. In the CNN-based scheme, the influence of neighbourhood pixels as input features is considered. However, a common practical problem based on deep learning classifiers is the class imbalance problem, which means that some neighbourhood classes have more examples in the training set than others. The overall accuracy measure is related to difficulties in the context of imbalanced data, and class imbalances can have an adverse effect on the training of deep learning classifiers [42]. This problem affects both the convergence of the training phase and the generalization ability of the model on the test set. In this study, due to the different neighbourhood sample sizes of the six categories, there is a certain class imbalance in the training and test samples in CNN compared with that in BPNN and RF, which may mislead the class label of training set and affect the CNN classification performance.
Although the three algorithm models have shown good performance in the classification of land cover changes in oceanic glaciers regions, the three algorithm models have different degrees of misclassification. The RF model incorrectly classifies "bare land" as "vegetation" and "vegetation" as "bare land", resulting in the phenomenon of unclear boundary between these classes of bare land and vegetation. Confusion between "vegetation" and "bare land" is also present in the CNN results, but it is not as strong as for the RF method. Both BPNN model and CNN model reduces the confusion between "ice" and "cloud" as well as "ice" and "shadow" compared to the RF method and improve the accuracy of these categories, but BPNN model still shows a tendency to misclassify "water". The artificial definition of thresholds between land cover categories for training samples may reduce the accuracy of land cover classification in oceanic glacial regions and thus, affects the analysis of land cover changes to a certain extent.

Driving Forces for Land Cover Change
Of the types of glaciers, oceanic glaciers are most sensitive to changes in temperature and precipitation, especially temperature changes [43]. An increase in temperature will cause glaciers to melt, and an increase in precipitation will increase the accumulation of glaciers. However, the increase in glaciers caused by the increase in annual precipitation is not enough to offset the increased ablation due to rising temperatures [44]. To understand the reasons for the changes in glaciers and vegetation, we downloaded the precipitation and temperature data of Bomê County from 1990, 2000, 2007 and 2016 from the China Meteorological Data Network (http://data.cma.cn/site/index.html) and calculated the total precipitation and average temperature one month before the image acquisition time in the study area (one month means 30 days), as shown in Figure 8. The precipitation in the study area showed a rapid increase from 1990 to 2000, with 95.6 mm of the monthly precipitation in 1990, increasing significantly to 148.0 mm in 2000. The monthly average temperature in the study area decreased from 10 degrees in 1990 to 8.78 degrees in 2000. Many studies have shown that an increase in temperature will cause glaciers to melt, and an increase in precipitation will increase the accumulation of glaciers [1,43]. Even if the temperature changes slightly, the glaciers may experience dramatic changes [11]. Consistently with the results of these scholars, the decrease in the average monthly temperature and the significant increase in the monthly precipitation can explain the dramatic increase in ice and decrease in the area coverage by vegetation from 1990 to 2000. After 2000, the average monthly temperature in the study area continued to increase at the same rate, reaching 9.91 degrees in 2007 and 11.93 degrees in 2016. For nearly three decades from 1990 to 2016, the average monthly temperature in the study area increased by approximately 2 degrees. In addition, the monthly average precipitation in the study area began to decrease, reaching 107.6 mm in 2016, but still higher than the original 95.6 mm in 1990. Chen et al. studied that the increase in glaciers caused by the increase in annual precipitation is not enough to offset the increased ablation due to rising temperatures [44]. Consistently with his results, glaciers decreased from 1990 to 2016 because the increase in glaciers caused by increased rainfall was less than that caused by temperature increases and confirmed temperature rise is the dominant factor in the rapid retreat of glaciers in this area [45]. From Figure 6d, the vegetation age began to gradually recover, and ice melted after 2000. The continuous decrease in precipitation and the relative increase in temperature can be attributed to this phenomenon.

Shortcomings and Prospects
Because glaciers in the southeast of Tibet are relatively least and the growth of vegetation can be effectively monitored in the summer, we generally prefer images acquired in summer to monitor land cover changes in southeast Tibet. However, due to the cloudy weather in summer, it is very difficult to obtain cloudless or low-cloud images in the summer in this area. Therefore, we selected four near-summer images from Bomê County (May 17, 1990, May 4, 2000, April 30, 2007, and May 24, 2016. Not only do they maintain a certain time span one year, but the time of image acquisition is spanning in one month. In future research, image fusion technology can be used to remove the cloud from the image to obtain a long-term continuous sequence of images to further analyze the land cover changes in the southeastern Tibet. Moreover, this study used Landsat images with a resolution of 30 m and a revisit period of 16 days. Due to the high cloud cover rate in the study area, Landsat images that can meet the cloud cover rate requirements are very limited, which results in coarse resolution. In future land cover classification in oceanic glacial regions, satellite data with high spatiotemporal resolution will obtain more information and perform more detailed analysis after data fusion. In addition, in terms of algorithms, although the CNN do not obtain the best classification results, in the field of deep learning, combining a CNN and other machine learning methods for land cover change analysis in oceanic glacial regions is also a very promising research direction [46].

Conclusions
The glaciers in southeast Tibet are the largest monsoon-type oceanic glaciers in China. As the temperature rises, the glaciers gradually melt, and the land cover in the glacial area also changes greatly, which has a significant impact on the stability of the ecosystem in Tibet. This study uses Landsat satellite images and three machine learning algorithms, namely, the RF, BPNN and CNN algorithms, to classify land cover in Bomê County. Three findings from the observations and algorithms are provided as follows.
First, the dominant land cover type was converted from ice to vegetation, and the landscape aggregation of ice decreased, indicating the influence of global warming. Second, in addition to the greatest transformation between vegetation and bare land, the probability of direct transformation between vegetation and ice increased when the elevation was higher than 4 km. Finally, although the RF and CNN algorithms have been widely applied in remote sensing image classification, this study shows the effectiveness of the BPNN algorithm in image classification for typical regions. Because clouds and mountain shadows influence the temporal and spatial quality of Landsat images for oceanic glacier mapping in Tibet, satellite data with high spatiotemporal resolutions are required in future studies, which should depend on the exploration of suitable data fusion algorithms.

Shortcomings and Prospects
Because glaciers in the southeast of Tibet are relatively least and the growth of vegetation can be effectively monitored in the summer, we generally prefer images acquired in summer to monitor land cover changes in southeast Tibet. However, due to the cloudy weather in summer, it is very difficult to obtain cloudless or low-cloud images in the summer in this area. Therefore, we selected four near-summer images from Bomê County (May 17, 1990, May 4, 2000, April 30, 2007, and May 24, 2016. Not only do they maintain a certain time span one year, but the time of image acquisition is spanning in one month. In future research, image fusion technology can be used to remove the cloud from the image to obtain a long-term continuous sequence of images to further analyze the land cover changes in the southeastern Tibet. Moreover, this study used Landsat images with a resolution of 30 m and a revisit period of 16 days. Due to the high cloud cover rate in the study area, Landsat images that can meet the cloud cover rate requirements are very limited, which results in coarse resolution. In future land cover classification in oceanic glacial regions, satellite data with high spatiotemporal resolution will obtain more information and perform more detailed analysis after data fusion. In addition, in terms of algorithms, although the CNN do not obtain the best classification results, in the field of deep learning, combining a CNN and other machine learning methods for land cover change analysis in oceanic glacial regions is also a very promising research direction [46].

Conclusions
The glaciers in southeast Tibet are the largest monsoon-type oceanic glaciers in China. As the temperature rises, the glaciers gradually melt, and the land cover in the glacial area also changes greatly, which has a significant impact on the stability of the ecosystem in Tibet. This study uses Landsat satellite images and three machine learning algorithms, namely, the RF, BPNN and CNN algorithms, to classify land cover in Bomê County. Three findings from the observations and algorithms are provided as follows.
First, the dominant land cover type was converted from ice to vegetation, and the landscape aggregation of ice decreased, indicating the influence of global warming. Second, in addition to the greatest transformation between vegetation and bare land, the probability of direct transformation between vegetation and ice increased when the elevation was higher than 4 km. Finally, although the RF and CNN algorithms have been widely applied in remote sensing image classification, this study shows the effectiveness of the BPNN algorithm in image classification for typical regions. Because clouds and mountain shadows influence the temporal and spatial quality of Landsat images for oceanic glacier mapping in Tibet, satellite data with high spatiotemporal resolutions are required in future studies, which should depend on the exploration of suitable data fusion algorithms.