Glacier Retreating Analysis on the Southeastern Tibetan Plateau via Multisource Remote Sensing Data

Accurate multitemporal glacier change investigations and analyses are lacking on the southeastern Tibetan Plateau (SETP). A combination of photogrammetry, optical remote sensing, and synthetic aperture radar datasets can accurately identify large-scale glaciers; In this article, glaciers in three periods on the SETP (1970s, 2000, 2020) were identified from multisource remote sensing data based on a deep learning method and manual visual interpretation, and multitemporal glacial inventory data from relatively high-frequency source imagery. Totals of 11648, 12993, and 11875 glaciers were identified in the 1970s, 2000, and 2020, with total areas of 13372.08 km2, 11692.31 km2, and 10612.94 km2, respectively. The general distribution of SETP glaciers was identified to be typical of alpine glaciers dominated by small-sized glaciers. The average elevation of glaciers was approximately 5000 m; the slopes were mostly lower than 40°, and the main aspect was southeast, followed by south and southwest. The glaciers retreated from the 1970s to 2020, and a total glacier area of approximately 2759.14 km2 was degraded during this time, with an average annual melting rate of 0.45% yr−1. Rising summer temperatures may be the driving force behind the continuous decline in the glacier area. Overall, the results obtained in this article showed relatively low uncertainty involved in the identification of glaciers compared to some previous studies. The results can provide accurate glacier information for glacier monitoring and modeling studies on the SETP.


I. INTRODUCTION
G LACIERS are sensitive to climate change, and the widespread shrinking of glaciers globally has been caused by climate warming in recent years [1], [2], [3]. The Qinghai-Tibet Plateau (QTP) is known as the third pole of the earth and the water tower of Asia, because its unique geographical conditions and abundant water sources [4], [5]. The QTP has Manuscript  also developed the largest mountain glaciers among the middleand low-latitude regions in the world [6], [7]. Throughout this region, the southeastern Tibet Plateau (SETP) has the most widely distributed marine glaciers, and many rivers, including the Yangtze River, Lancang River, and Brahmaputra, originate from the QTP and flow through the SETP [8]. The study of glacier changes is of great significance to hydrology and water resource research and to regional climate change investigations. Glaciers on the SETP have experienced significant mass losses in recent decades, and their mass loss rate was approximately three times larger than that of the entire QTP [1], [9], [10]. At the same time, similar to the continuous retreat of the entire QTP glacier, glacier retreat in this region has also been reported by some researchers [5], [11].
Regarding the research on glaciers on the SETP, few studies have investigated the glacier area changes by focusing only on some typical glaciers in the region [12], [13], [14]. Research on the mass balance of glaciers, especially large-scale glaciers, is relatively mature, combining laser altimetry (ICESat/ICESat-2, Cryosat/Cryosat-2), digital elevation model (DEM) and gravity satellite data such as GRACE and GRACE follow-on data [1], [9], [15]. However, the glacier inventory information that is often used for large-scale mass balance investigations is derived from satellite imagery of varying quality acquired in different time periods, such as Randolph Glacier Inventory (RGI) data, which involve uncertainty. Therefore, obtaining multitemporal glacial inventory data from relatively high-frequency source imagery may allow researchers to cope with such uncertainty.
For the identification of glaciers on the SETP, the first Chinese Glacier Inventory (CGI-1) started in 1978, and the inventory information was completed in 2002 [16]. The glacier outlines were sketched based on aerial photographs and digital topographic maps obtained from the 1950s to the 1980s. However, due to this long temporal span, the inventory cannot represent the current status of glaciers. In the second Chinese Glacier Inventory (CGI-2), although Landsat-5 TM images were used as the main data source from 2004 to 2011, it was difficult to obtain high-quality optical satellite images since the inventory information for the SETP follows CGI-1 [17]. In addition, global glacier inventory information, including information characterizing the SETP region, such as the China-region glaciers in the RGI, has been digitized from scanned digital copies of 80 glacier inventory maps produced during the CGI compilation process [6], [8].
Large-area glaciers can be identified from cloud-and snowfree high-resolution optical remote sensing data. The Landsat This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ and Sentinel-2 satellites can be used as the main data sources, and glaciers can be easily identified by artificial visual interpretations or the threshold method [13], [18], [19]. However, the SETP is often affected by the South and East Asian monsoons [20]; Summer precipitation accounts for more than half of the annual precipitation, and the large-scale snow cover in winter and cloud cover in summer limit the application of optical remote sensing data over the SETP [21]. The systematic and large-scale inventory of glaciers on the SETP was launched mainly by the RGI and the first and second glacier inventories in China, while the second glacier inventory was performed based on the first glacier inventory due to the lack of available optical images of the SETP [6], [8], [17]. Synthetic aperture radar (SAR), which is not affected by weather conditions, can supply optical remote sensing data; thus, SAR as the data source has been used to identify glacier in recent years [22], [23], [24].
In view of the long-time span of the images used to construct large-scale glacier inventories of the study region, the multitemporal glacier inventory on the SETP are lacking. Combining multiple optical remote sensing and SAR images makes it is possible to obtain glacier inventories over multiple periods. For optical imagery, the earliest Landsat satellite was launched in the 1970s. Currently, Landsat 1/2/3/5/7/8/9 are freely available and as the main data source are used to identify glaciers [6], [18], [25]. Sentinel-2 can obtain images with a resolution of up to 10 m since 2017. Compared to Landsat's 30-m-resolution data, Sentinel-2 has a great advantage; therefore, Sentinel-2 is often used as the main data source for glacier identifications [19], [26]. SAR data were developed at the end of the 20th century. The main SAR data sources commonly used for glacier identifications include the ERS (ERS-1, ERS-2) [27], ENVISAT ASAR [28], ALOS PALSAR [29], TerraSAR [30], and Sentinel-1 [31] satellites. In studies on large-scale glaciers in the early period of related research, such as before 1980, the first-generation optical Landsat multispectral scanner system (MSS) image data with a 60-m resolution could be used. However, for high-resolution analyses, public digital aerial photography images, such as keyhole images (KH-4, KH-9) with a resolution of 6-9 m, may be used for classification tasks [13], [32].
The rapid development of Convolutional Neural Networks (CNNs) enables it to be widely used in image classification [33], [34]. Since then, the emergence of Visual Geometry Group (VGG), GoogLeNet, and residual neural network (ResNet) have greatly improved the classification accuracy of images. In order to obtain deeper semantic information, more advanced networks have been developed, such as Unet, SegNet, DeepLab. CNNs is widely used to glacier identification due to its excellent texture feature extraction ability [35], [36]. Several studies have shown that CNNs have higher accuracy in glacier identification [35], [37]. There are relatively few applications of deep learning in the study of mountain glaciers and optical remote sensing is easily affected by cloud. Moreover, it is not conducive to the identification of debris-covered glaciers because the similar spectral characteristics of debris-covered glaciers and surrounding rocks, so it is still challenging to extract glacier boundaries quickly and semiautomatically based on remote sensing images [38]. Deep learning can automatically extract features of different depths from the input feature image through CNNs and combines the deep features obtained by semantic segmentation with the interference and planned shallow features of the glacier surface to provide automatic, fast and high-precision glacier identification results [37], [39].
To obtain an inventory information of glaciers for more than two points at each time step on the SETP and analyze the results, the following step-by-step methodology was proposed and applied in this article.
1) For the 1970s, KH-9 images were used as the main data source, and Landsat MSS images were used as the supplementary data source to obtain a glacier inventory by artificial visual interpretation combined with the outlines of the RGI 6.0 and CGI1.0 datasets. 2) For 2000, Landsat TM/ETM+ images were preferentially used, and ERS-2 images were used in cloud-covered areas.
The preliminary glacier boundaries were sketched by a deep learning method, and the glacier inventory in 2000 was obtained by adding artificial corrections based on the true images. 3) To obtain the inventory for 2020, the same method as that used for 2000 was followed, but the Sentinel-2 and Landsat OLI images were used as the main data sources, while Sentinel-1 images were used for the cloud-covered areas. 4) Glacial changes were then analyzed for the period from the 1970s to 2020 based on the three-period glacier inventory. 5) The uncertainty involved in the determined glacier area was assessed using the buffer method and empirical formula, combined with a deep learning-based classification method. 6) Finally, considering the air temperature and precipitation factors, the possible drivers of glacier changes were explored and investigated.

II. STUDY REGIONS
The SETP is located on the southeastern Qinghai-Tibet Plateau (QTP) (see Fig. 1), with an average elevation of more than 4000 m. The development of glaciers on the SETP is dominated by marine glaciers, and most of the glaciers are distributed in the western Hengduan Shan, eastern Nyainqentanglha and eastern Himalaya. The RGI 6.0 dataset contains 11170.94 km 2 of glaciers on the SETP, accounting for approximately 1/5 of all QTP glaciers [6], [8]. In recent years, the overall glacier mass has been faced with accelerated loss, and the amount of mass lost here has been higher than that on the QTP [40]; such a mass loss has led to an increase in the number of glacial lakes on the SETP [41], [42].
The SETP is affected by the southeast monsoon, the South Asian monsoon and the westerlies. The precipitation in this region is generally higher than that in other QTP regions. In the southeastern region with lower elevations, the annual precipitation can exceed 1000 mm [40]. The observed data from seven automatic meteorological stations in the study area since the 1970s were gathered from the Meteorological Information and Network Center of the Tibet Meteorological Bureau for use in this article. These stations are concentrated in  the Nyainqentanglha region, namely, Chamdo, Lhorong, Pome, Paksho, Nyingtri, Zogong, and Zayul. According to the monthly average temperature and precipitation records collected from 1970 to 2020, precipitation on the SETP is concentrated in summer, with only a small amount of precipitation occurring in winter (approximately 5% of the annual precipitation, Fig. 2).

A. KH-9 Images
The KH-9 satellite (the Hexagon mission) was a spy reconnaissance satellite launched by the United States in the 1970s. The KH-9 data are panchromatic images with resolutions of 6-9 m [43]. Sixty-two KH-9 images (Table S1) taken from 1973 to 1976 were used to identify glaciers in the 1970s as the main source of the data. The data were ordered and downloaded from USGS Earth Explorer (https://earthexplorer.usgs.gov/, last access: August 23, 2022).  (Table S3) and 26 Landsat-7 images (Table S4) and 3 bands (SWIR, NIR, and Red) with a resolution of 30 m were used for the 2020 glacier identification. From Landsat-8, which was launched on February 11, 2013, 21 images (Table S5) with NIR, Red, and Green bands and a spatial resolution of 30 m were employed as supplementary data to identify the glacier boundaries over some areas in 2020. The panchromatic band 8, with a resolution of 15 m, was also used as a reference orthoimage for the absolute orientation of the KH-9 image. All Landsat images could be obtained from the USGS Earth Explorer (https://earthexplorer.usgs.gov/, last access: August 23, 2022). 2) Sentinel-2: Sentinel-2 is a high-resolution multispectral imaging satellite carrying the Multispectral Imager. The Sentinel-2 satellite is 786 km high and has a 290-kmwide swath, with 13 spectral bands ranging from visible light and near-infrared to shortwave infrared wavelengths, with a maximum resolution of 10 m [19]. A total of 50 Sentinel-2 images (Table S6) taken from July to September 2020 were downloaded from Copernicus (https:// scihub.copernicus.eu/dhus/#/home, last access: August 23, 2022) and used as the main dataset for the 2020 glacier identification.

C. SAR Satellite Images
Optical images are easily affected and limited by clouds, especially on the SETP and more specifically in the rainy summer season. However, SAR images are not affected by clouds and can be used in areas where glaciers cannot be identified in optical images.
1) ERS-2: The ERS-2 satellite is the second resource remote sensing satellite of the European Space Agency (ESA) and was launched in 1995. The spatial resolution of the ERS-2 SAR data is 26 m, the image width is approximately 100 km, and the revisit period is 35 days. Thirtysix ERS-2 SAR images (Table S7) obtained from the European Space Agency (https://esar-ds.eo.esa.int/oads/ access/collection, last access: August 23, 2022) could be used to compensate for the missing parts of the Landsat-5 and Landsat-7 optical images to obtain the glacial boundaries in 2000. 2) Sentinel-1: Sentinel-1 is a near-polar sun-synchronous orbit satellite launched by the ESA in April 2014. The Sentinel-1 mission provides data from a dual-polarization C-band SAR instrument at 5.405 GHz (C band) and has four imaging modes: stop map, interferometric wide swath, extra wide swath, and wave, as well as both up and down orbit data. The satellite has a 12-day revisit period and provides all-day, all-weather radar imagery for land and marine services [16]. Forty Sentinel-1 images (Table  S8) were used in this article to identify the glaciers in 2020 from ESA (last access: August 23, 2022).

D. Previous Glacier Inventories
The RGI is a complete inventory of global glaciers published by Global Land Ice Measurements from Space. The RGI glacier inventory has been released in 6 versions. The first version was released in February 2012, and the latest version was released in July 2017; the latest version includes the second glacier inventory data in China. In SETP, RGI 6.0 includes 12 166 glaciers. Optical images from 1999 to 2013 were used for glacier identification: 5652 glaciers from 1999-2004 and 6514 glaciers from 2005-2013 were counted [6]. RGI 6.0 is usually used as a reference guideline for the glacier identification process and is freely available via http://www.glims.org/RGI (last access: August 23, 2022).
The CGI-1 is a part of the Glacier Inventory Project conducted by the Institute of Environment and Engineering in Cold and Arid Regions of the Chinese Academy of Sciences from 1987 to 2002. The dataset covers all of China and was published by the China Glacier Information System in 2004. The CGI-1 inventory was constructed based on aerial topographic maps and aerial photographs taken from the 1950s to the 1980s through manual interpretation. On the SETP, CGI-1 includes 12 197 glaciers [16], [17]. The CGI-1 data can be obtained from the National Glacier and Permafrost Desert Science Data Center (http://www. ncdc.ac.cn/portal/, last access: July 4, 2022) and can be used as auxiliary data in the glacier identification process.

E. DEM and Climate Data
The NASA DEM (NASA DEM) with a 30-m resolution was used to extract topographic information on the glacier surfaces via https://www.earthdata.nasa.gov (last access: August 23, 2022), including the slope, aspect and elevation information of each glacier. In addition, a terrain correction was applied to the ERS-2 and Sentinel-1 data and an absolute orientation step was applied to the KH-9 data to facilitate the use of the NASA DEM. The NASA DEM represents a significant improvement among the available three-arc second Shuttle Radar Topography Mission DEM, and since it integrates the GDEM with other data sources, the coverage and accuracy of the resulting data are greatly improved.
Climate factors such as temperature and precipitation were also used to analyze the driving factors of glacier area and mass changes. Monthly air temperature (recorded at 2 m above the ground surface) and precipitation time series from the European Center for Meteorology and the ERA5-Land forecast reanalysis product were extracted for the study area from 1973 to 2020. The dataset was downloaded from https://cds.climate.copernicus.eu/ (last access: August 23, 2022) and regridded to a regular latitudelongitude grid of 0.1 × 0.1°.

A. Glacier Inventory in the 1970s
Due to the missing parameters and serious distortion of the KH-9 data, the downloaded unprocessed images could not be conformed to the actual geo-reference, and the conventional control point correction method could not achieve large-scale geo-referencing results [43]. Therefore, existing mature remote sensing software cannot be used directly to register or generate DEMs. KH-9-based DEMs and ortho-images can be obtained using the open-source software Micmac [44]. The operation process can be divided into data preprocessing, relative orientation, and absolute orientation (see Fig. 3). These operations were carried out in the Python 3.8 environment. The automatic extraction of ground control points was performed using Landsat-8 ortho-images and 30-m NASA DEM data to obtain the reference coordinates and elevation values in the absolute orientation.
Distortion-corrected and geocoded KH-9 ortho-images were used as the primary data sources to extract glacier outlines in the 1970s. Since the KH-9 images were grayscale images, conventional threshold classification methods could not be used for automatic extraction. Moreover, it was difficult to distinguish the surface moraine area from the surrounding rocks, which further increases the difficulty of extracting the glacier boundaries. Therefore, the manual visual interpretation method was used for glacier identification. To improve the visual interpretation accuracy, the 5% linear stretching method [41] was used to enhance the KH-9 images of the glaciers and surrounding bare rocks [see Fig. 4(a)]. This glacier outline identification process was performed in reference to RGI 6.0 and CGI 1.0.
Most glaciers were identifiable from the KH-9 images, but since some of the KH-9 images were acquired during winter when some areas were covered by clouds, the Landsat MSS images were used as a secondary data source to identify glacier outlines using manual visual interpretation [see Fig. 4(b)].

B. Glacier Inventories in 2000 and 2020
To identify glaciers in 2000, Landsat-5 data were preferentially used as the main data source, and Landsat-7 data were used as supplements for the areas that were not covered by Landsat-5. ERS-2 SAR data were used as surrogate data for areas that were not covered by the optical images. For the 2020 glacier identification, Sentinel-2 images were used as primary data, accompanied by Landsat-8 data and Sentinel-1 SAR data for the areas that were not covered by the optical imagery.
The optical image processing steps included cloud removal, band composition (red, green and blue) and the calculation of the normalized difference snow index (NDSI). However, for SAR images, the data processing methods contained the following steps: performing orbit correction on the original image, performing radiometric calibration on the corrected image, and converting the backscatter of the image into absolute radiance. Finally, image filtering should be performed, and the filtered image should be georeferenced to eliminate the geometric distortion. In this article, the NASA DEM was used to correct the ERS-2 and Sentinel-1 images.
A deep learning method was used to automatically identify glacier outlines. The UNET network has a U-shaped structure, mainly including an encoder and a decoder [45]. The encoder reduces the dimension of the image and extracts features through convolution and pooling. The combination of layer features and deep features is more conducive to glacier extraction [37]. The VGG16 convolutional neural network model consists of 13 convolutional layers and 3 fully connected layers. In this article, VGG16 is used as the encoding part to extract glacier features, so only the first 13 convolutional layers are selected, and the convolution kernel size is 3 × 3. The stride is 1, filled in the same way, and the activation function is the relu function. VGG uses smaller convolution kernels and deeper convolution layers to replace neural networks with shallower layers and larger convolution kernels, and can extract deeper features improve the recognition ability of the network [39]. The 13-layer convolution can be divided into 5 types of convolutions, the number of layers of each type of convolution is 2, 2, 3, 3, 3, and the number of convolution kernels of each type is 64, 128, 256, 512, 512, and the number of layers and the maximum pooling layer is used to connect the layers, and its specific structure is shown in Fig. 5.
For optical remote sensing data, including Landsat 5/8 and sentinel-2, the parameter information of the three bands (SWIR, NIR, and Red) is used as the sample part of the deep learning network, and then the histogram equalization process is performed to make the image information richer, the features are more obvious. For the classification sample selection of SAR data, in 2000, two input parameters are used as the sample part of the deep learning network includes the backscatter under the ERS-2 VV polarization and the DEM. Similarly, in 2020, three input parameters includes the VH and VV information of Sentinle-1A and the DEM are used as the sample part of the deep learning network (see Fig. 6). Based on the RGB composition of optical images, NDSI and RGI 6.0, visual interpretation was used to identify 2000 glaciers in the study area, and the 2000 glaciers as training datasets in deep learning; Select the 2000 glacier boundaries for binary segmentation as the label part of the network. The labels of the glacier include pure glacier and debris-covered glacier.
The size of the sample is cut to 512 × 512, and there are a total of 5000 samples. The samples and labels are divided into training datasets, verification datasets, and test datasets according to 6:2:2, and data enhancement processing is performed on the training datasets and verification datasets. By expanding the dataset to enhance the generalization of the model. Based on the deep learning semantic segmentation network VGG16-UNET, the training datasets and the verification datasets are trained, and finally applied to the test datasets to obtain the binary classification results of glaciers and nonglaciers (see Fig. 7).

C. Manual Correction Principle
To obtain more accurate glacier outlines, we manually corrected the three-stage glaciers. The glacier identification method for the 1970s was based on visual interpretation. Hence, we manually corrected the preliminary glacier outlines obtained by the deep learning method in 2000 and 2020. The following unified manual correction standards were applied. For the visual interpretation of glaciers in the 1970s, the outline of each glacier was modified on the basis of RGI 6.0 and CGI 1.0, and the KH-9 image was preferentially used. If the KH-9 image was covered by a large amount of snow or there was no KH-9 image of the area, a correction based on the Landsat MSS image was applied. If neither the KH-9 or Landsat MSS image could not be judged, the RGI 6.0 or CGI 1.0 information would be retained. For 2000, based on RGI 6.0 and the true-color composition of the Landsat-5 and Landsat-7 imagery, if the identified glacier area was covered by the RGI, it would be corrected; if there was no RGI coverage but the identified glacier area was greater in area than 0.1 km 2 , it was also corrected. If both the Landsat-5 and Landsat-7 imagery in the identification area were covered by clouds, the outlines were retained. In addition, if the glacier area was less than 0.1 km 2 , it was deleted. For 2020, the applied steps were the same as those used for the manual correction in 2000, but the optical images were based on the Sentinel-2 or Landsat-8 imagery.
In the process of identifying glaciers by VGG16-UNET, the training datasets includes pure glaciers and debris-covered glaciers. Therefore, the results only include glacier areas and nonglacial areas, and the glacier area includes the debris-covered glacier. Since there is no available optical image for the glacier area identified by SAR, the identification result of deep learning is retained without any correction. For glacier areas identified by optical images, NDSI and Google Earth images are used to manual correction. The main part of the glacier is corrected according to the true color image, and the glacier with debriscovered is corrected with NDSI. The correction process mainly refers to RGI 6.0. For 2000 and 2020, if the target glacier is partially covered by clouds on the optical image, then refer to the NDSI image and Google Earth current and historical images to correct it (see Fig. 8). Only when optical images, NDSI and Google Earth are not available, no corrections are made.

D. Glacier Division and Inventory Information
The identified glaciers were divided based on the dividing criteria for the RGI; then, the number of glaciers could be counted. In addition, in the 1970s, for the glacier area not covered by the RGI, the division principle of the CGI was adopted. On the other hand, for 2000 and 2020, the glaciers that were identified in this article, but not in the RGI, were not divided.
After determining the outline of each glacier, the outlines of glaciers in the 1970s, 2000, and 2020 were stored in Shapefile format, and the ID of each glacier was unified based on the latitude and longitude. The topographic information of each glacier, including its average slope and aspect (both in degrees) and its maximum, minimum, and mean elevation, was calculated from the DEM. The DEM in the 1970s was from the DEM generated from KH-9 imagery, and the ungenerated areas in the 1970s, 2000, and 2020 refer to the NASA DEM. The glacial area (km 2 ) and perimeter (m) were also calculated using the equal-area projection tool (WGS84 coordinate system, Asia North Albers Equal Area Conic) in ArcGIS. Further information was added to each glacier in the attribute table as follows: the name of mountain system, the satellite used, the acquisition time of the satellite, the source of the glacier elevation, the acquisition time, and the glacier name (from the RGI 6.0).

E. Uncertainty Assessment
In remote sensing-based identifications of glacier area, uncertainty arises from the characteristics of the source data and the errors caused by the later manual correction. The former includes cloud cover, seasonal snow cover, shadows, and glacier debris, all of which can lead to misclassifications of remote sensing images when using different kinds of identification methods [18].
For optical image identification, we superimposed the reported results by other studies and added a few manual interpretations and inspections in the later stage. Therefore, the source of error was mainly caused by the image resolution. The buffer method was used to evaluate the uncertainty due to the resolution [18], using the glacier perimeter multiplied by 1/2 the image resolution as follows: where Uncertainty is the error area caused by the resolution, L is the perimeter of the glacier, and R is the pixel resolution. In addition, we employed the glacier error model used by the RGI 6.0 to determine the area error for each glacier as follows [6]: where e(s) represents the error under the glacier area of s, e 1 is the fractional error under 1 km 2 , which is 0.039, p is the correction factor, which was set to 0.7 in this article, and k = 1.
For the glacier classification of SAR data by deep learning, the source of uncertainty also depends on the uncertainty caused by the deep learning method. The overall accuracy of the modeling results was used to evaluate the accuracy of the model as follows where TP indicates a true identification of a positive sample, TN indicates a true identification of a negative sample, FP indicates a false identification of a negative sample as a positive sample, and FN indicates a false identification of a positive sample as a negative sample.

A. Glacial Characteristics in the 1970s, 2000, and 2020
In the 1970s, 2000, and 2020, totals of 11648, 12993, and 11875 glaciers larger than 0.01 km 2 were identified on the SETP, with areas of 13372.08 km 2 , 11692.31 km 2 , and 10612.94 km 2 , respectively (see Table I). The number and area of glaciers in the three subregions were counted separately. In the 1970s, the numbers of glaciers in the Hengduan Shan, Nyainqentanglha and Eastern Himalaya (the same below) were 2034, 7191, and 2423, respectively, with areas of 1388.33 km 2 , 8943.98 km 2 , and 3039.77 km 2 , respectively. In 2000, 2072, 8165, and 2756 glaciers were identified, with areas of 1355.07 km 2 , 7407.28 km 2 , and 2929.96 km 2 , respectively. In 2020, 1981, 7295, and 2599 glaciers were identified, with areas of 1238.41 km 2 , 6657.58 km 2 , and 2716.95 km 2 , respectively. More than 60% of glaciers were distributed in Nyainqentanglha, and Hengduan Shan included the fewest number of glaciers and the smallest glacier area (<15%).
The numbers and areas (see Table II) of glaciers of different sizes were further counted. In general, the number of small-area glaciers (<0.5 km 2 ) on the SETP dominated the overall number of glaciers, but their area was relatively small. In the 1970s, 59.23% (76.54%) of glaciers with areas less than 0.5 km 2 (1 km 2 ) accounted for only 0.76% (11.04%) of the glacier area. In 2000 and 2020, glaciers with areas less than 0.5 km 2 (1 km 2 ) accounted for more than 69% (83%) of the total number of glaciers but accounted for only approximately 1.6% (14%) of the glacier area. Only 3.76%, 2.88% and 2.82% of glaciers larger than 5 km 2 were identified in the 1970s, 2000, and 2020 (the same below), respectively, but accounted for 42.92%, 43.72%, and 42.82% of the total glacier area, respectively. Therefore, glaciers with areas between 1 and 5 km 2 accounted for 19.70%, 13.72%, and 13.77% of the total number of glaciers in the 1970s, 2000, and 2020, respectively, and for 35.43%, 31.25%, and 31.68% of the total glacier areas in these years, respectively. This difference in number and area is typical for alpine glaciers [19].
The terrain characteristics, including the average elevation, average slope, and average aspect, of each glacier were calculated (see Fig. 9). In general, glaciers have a high elevation proportion at approximately 5000 m, and the highest proportion of slopes are between 10°and 20°. This distribution shows that the glaciers on the SETP face southward, with the highest proportion facing southeast, followed by south and southwest.

B. Glacial Changes
From the 1970s to 2020, the total area of the glaciers on the SETP decreased by 2759.14 km 2 (−20.63%), with a decrease   Table III). Since most of the KH-9 data used in the 1970s were from 1973-1975, when calculating the areal changes, the areal changes of the 1970s were based on the year 1974. Among them, the glacier area decreased by −1679.77 km 2 (−12.56%) and by −1079.37 km 2 (−9.23%) from the 1970s to 2000 and from 2000 to 2020, respectively, and the reduction rates were −64.61 km 2 yr −1 (0.48% yr −1 ) and −53.97 km 2 yr −1 (0.46% yr −1 ), respectively. Glaciers retreated faster in the former period than in the latter, but the difference was not significant. Among the three regions, Nyainqentanglha contributed the largest glacier retreating area and rate, which was higher than the ablation rate of the overall region. While in Hengduan Shan and Eastern Himalaya, the ablation rates were lower than that of the whole study area, at 0.24% yr −1 and 0.23% yr −1 , respectively, the ablation rates from 2000-2020 (0.43% yr −1 and 0.36% yr −1 ) were greater than those from 1970-2000 (0.09% yr −1 and 0.14% yr −1 ). According to Fig. 10 and the characteristics of the glaciers presented in Tables I and II, it can be seen that in different periods, the area and number of glaciers larger than 0.5 km 2 continuously decreased from the 1970s to 2020, while glaciers smaller than 0.5 km 2 have shown an increasing trend first followed by a decrease, showing an overall increase, and the number and area of these glaciers have increased by approximately 10% and 3%, respectively. The changes in the areas of glaciers of different sizes in different periods and their contributions to glacial retreat were calculated and are presented in Fig. 10.
The area of glaciers larger than 10 km 2 continued to decrease, the area of glaciers from 2000 to 2020 decreased more than that from the 1970s to 2000, and the areal reduction of these glaciers accounted for approximately 30% of the overall glacial area reduction. In the glacial retreat process, glaciers experience areal reductions, and the glaciers themselves are further split until they disappear. In general, the increased number of glaciers from the 1970s to 2000 could be linked to the increase in the number of glaciers with areas <0.5 km 2 , and the splitting of glaciers into two or more glaciers could be the main reason for this phenomenon. From 2000 to 2020, both the area and the number of glaciers decreased to varying degrees, and the glaciers that split from the 1970s-2000 retreated further or even disappeared from 2000-2020. Fig. 11 shows an example of glacial retreat. The glacier area is decreasing, and glacier retreat is observed at the glacier toe. Smaller glaciers could disappear entirely under warming temperatures.

A. Comparisons With Previous Studies
Our findings provide information on the distribution of glaciers in southeast Tibet in the 1970s, 2000, and 2020, compared with previous studies, our study uses multiple remote sensing data sources of the same period to reduce the uncertainty caused by changes in glaciers. In 2020, for the first time, higher-resolution sentine-1 and sentinel-2 images were used to identify glaciers in SETP.
Most of the aerial images used in the first glacier inventory in the SETP were taken from 1956 to 1986, with a time span of up to 30 years, and the glacier area in the original CGI-1 was based entirely on manual measurements obtained with a planimeter, so these data were subject to several false reports and the shadow influence. The glacier area on the SETP was overestimated by the original CGI-1 [16], [17]. The CGI-1 identified 12 198 glaciers in southeastern Tibet with an area of 17 120.53 km 2 ; this area is much larger than the result obtained in this article based on KH-9 images from 1973 to 1976 (see Table IV).
From the latest version of the global glacier inventory product RGI 6.0, 12166 glaciers could be identified with an overall area of 11 170.94 km 2 on the SETP [6]. Compared to the glacier area obtained for 2000, the area has decreased by 521.37 km 2 . The main reason is that the time span of the RGI 6.0 image source was large, spanning from 1999 to 2013, while most of the images from Landsat-5, Landsat-7, and ERS-2 used in this article were from 2000. Our results were also plausible as the glaciers continue to decline.
The areal reduction rate of glaciers on the SETP was 0.45% yr −1 since the 1970s, showing agreement with the results  provided in some previous studies for the SETP. For example, Liu et al. [41] stated that the areas of 141 glaciers in the Bugyai kangri area in the Nyainqentanglha region decreased by 30.44 km 2 from 1981 to 2013, with a reduction rate of 0.48% yr −1 . Wang et al. [14] identified a total glacier area Although an ablation rate of 0.24% yr −1 for Hengduanshan was obtained in this article, the loss rate was as high as 0.43% yr −1 between 2000 and 2020, compared to only 0.09% yr −1 before 2000. Pan et al. [12] analyzed the glaciers in the Gongga Mountains area under the Hengduan Mountains and obtained a glacier ablation rate of 0.26% yr −1 from 1966 to 2009, which was close to the overall results obtained for the Hengduan Mountains in the current study.

B. Uncertainty in Glacier Identification
To obtain the glacier outlines from the optical images, a large number of manual corrections were applied; therefore, uncertainties could arise from the resolution parameters of the remote sensing images used and the errors in the sample selection process. To identify SAR images by deep learning, the uncertainty could be linked to the accuracy of the classification model and the errors caused by the resolution of the remote sensing images when smoothing the glacier boundaries [18].
The overall accuracies of the VGG16-UNET when classifying glaciers from ERS-2 and Sentinel-1 imagery were 82% and 93%, respectively. The results were comparable to the accuracy of some other deep learning methods used for glacier identification [31], [36], [37], [38]. In 2000, limited by the quality of the data themselves, the only data input to the deep learning network was VV-band data from ERS-2, while the Sentinel-1 data used in 2020 comprised two input bands, VV and VH. The combination of multiple bands can often lead to better results than when only one band is used to identify and extract glacier features.
To quantify the uncertainty linked to the perimeters of the glaciers, the resolution of the source images, and the accuracy of deep learning for SAR identification, the buffer area method and error model [18], [19] were applied to the glaciers identified in the 1970s, 2000, and 2020 (see Fig. 13).
The area errors caused by the image resolution were 364.13 km 2 , 756.14 km 2 , and 322.63 km 2 in the 1970s, 2000, and 2020, respectively, accounting for 2.72%, 6.47%, and 3.03% of the corresponding glacier areas (see Table V). According to the uncertainty analysis, in the 1970s and 2020, the error caused by the resolution was lower than that reported by some previous studies on glacier inventories [6], [13], [17], mostly due to the high-resolution images used in this article. Most of the identification results in the 1970s came from KH-9 data. In 2020, nearly 70% of the glaciers were interpreted from Sentinel-1 and Sentinel-2 data. The image sources of some glacier inventory products mainly include Landsat imagery with a resolution of 30 m. In 2000, because the resolution of the images used was lower than those used in the other two periods, the error was greater than that in the 1970s and 2020, but the proportion was comparable to those reported in other studies [6], [13], [19]. In addition, the errors of the calculated areas extracted by the glacier error model were 364.13 km 2 , 756.14 km 2 , and 322.63 km 2 in the 1970s, 2000, and 2020, respectively, accounting for 2.92%, 3.04%, and 3.05% of the glacier area, and the accuracy was close to that reported in RGI 6.0 for China. Among them, in the 1970s and 2020 results, the errors were close to the error areas caused by the buffer method. In 2000, apart from the effect of low-resolution data, due to an increase in the number of glaciers, the peripheral area of the buffer zone increased, which could have increased the uncertainty of the buffer method.

C. Climate Factors Affecting Glacier Changes
Changes in temperature and precipitation are usually considered to be the main climatic drivers of glacial retreat [5], [7], [13]. Since most of the glaciers on the SETP are typical marine and summer-accumulating glaciers, the climatic factors from May to September can best explain the glacial changes [14], [41]. Air temperature and precipitation data from the ERA-5 Land dataset were used to extract climate factor changes on the SETP. The average annual temperature in the Nyainqentanglha region was much lower than that in the Hengduan Mountains and the Eastern Himalaya [see Fig. 14(a)], so the numbers and areas of glaciers were much larger in these regions than in the other two regions, and the rate of glacial ablation was slow. The annual average temperature in southeastern Tibet showed a significant warming trend from 1973 to 2020. Among them, the heating rate of the Hengduan Mountains was significantly higher than that of Nyainqentanglha and the Eastern Himalayas, while the warming rates of Nyainqentanglha and the Eastern Himalayas were almost the same. The annual precipitation has shown a slight downward trend since 1973 [see Fig. 14(b)]. Throughout the study area, the Eastern Himalayan region was most affected by the South Asian monsoon due to its proximity to the southern area, and its annual precipitation was higher than those of the other two regions.
An increase in temperature causes the surface of a glacier to absorb more heat radiation, thus reducing the albedo of the glacier and accelerating glacial melting and mass loss [5], [46]. An increase in temperature may be the main factor driving the retreat and ablation of glaciers on the SETP.

VII. CONCLUSION
Multisource remote sensing data, visual interpretation and a deep learning method were used in this article to identify glaciers on the SETP in the 1970s, 2000, and 2020. The overall uncertainty was quantified by combining the uncertainties linked to the image resolution and by using a deep learning classification method. The characteristics of the glacier area changes from the 1970s to 2020 were obtained, and the driving forces of climatic factors on glacier changes were discussed. The following main conclusions can be drawn based on the results obtained in this article.
1) A total of 11 648 glaciers with a total area of 13 372.08 km 2 were identified using the KH-9 and Landsat 1-3 datasets in the 1970s. A total of 12 993 glaciers with a total area of 11 692.31 km 2 were identified from ERS-2, Landsat-5, and Landsat-7 data in 2000. A total of 11 875 glaciers with a total area of 10 612.94 km 2 were identified from the Sentinel-1, Sentinel-2, and Landsat-8 data in 2020.
2) The distribution characteristics of the glaciers were typical of mountain glaciers, with a large number of small-sized glaciers but a low areal proportion of these small glaciers. The average elevation of the glaciers was close to 5000 m, the slopes were mostly lower than 40°, and most glaciers faced southeast, south, or southwest. Nyainqentanglha showed the largest glaciers in the study area, followed by the Eastern Himalaya and Hengduan Shan.
3) The area of glaciers continued to decrease from the 1970s to 2020, with a total degraded glacier area of 2759.14 km 2 , corresponding to an average annual loss of 59.98 km 2 and a melting rate of 0.45% yr −1 . The increase in the number of glaciers from the 1970s to 2000 was mainly caused by the splitting of glaciers. The decreases in the number and area of glaciers from 2000 to 2020 were due to the continuous decline in glaciers and the disappearance of some small glaciers. 4) Due to the application of improved and high-resolution images for the 1970s and 2020 and the reliable deep learning classification method, the uncertainty of the glacier area estimation results was less than that reported in some previous studies on large-scale glaciers. 5) The air temperature on the SETP has shown a significant warming trend since 1973, which may be the main driving force for glacial retreat in this region.