Monitoring Rock Desert Formation Caused by Ice–Snow Melting in the Qinghai-Tibet Plateau Using an Optimized Remote Sensing Technique: A Case Study of Yushu Prefecture

Abstract

The rapid and effective identification of the spatial distribution of rock deserts in ice–snow melting areas can provide useful information for the prevention of natural disasters. In this study, Landsat TM/OLI were used to track the evolution of snow lines in the snow melting areas of Yushu Tibetan Autonomous Prefecture of the Qinghai-Tibet Plateau (hereinafter referred to as Yushu Prefecture) during the last 30 years. A total of seven extraction schemes were used to classify the rock desert in the ice–snow melting area by remote sensing. Our results show that: (1) The accuracy of the multi-index factor compound analysis and object-oriented classification compound method was the highest. This method can provide a rapid and efficient reference scheme for rock desert extraction in the plateau ice–snow melting area. (2) The combinations of two single methods with higher accuracy can further improve the total accuracy. If a single method with lower accuracy is involved in the multi-method fusion, the accuracy of the method with lower accuracy can be improved. (3) In the past 30 years, there has been a large amount of ice–snow melting in Yushu Prefecture. The ice–snow melting area accounts for 53.78% (1451.04 km²) of the ice–snow area, and the bare rock and bare stone in the rock desert account for 63.77% of the total area of the ice–snow melting area. Bare sand and bare soil area account for 30.27% of the total area of ice and snow melting area.

1. Introduction

With the rise in global temperature, the last 30 years have been the warmest 30 years in the past 1400 years [1]. High latitude and high-altitude areas are obviously affected by global climate change [2]. Since the 20th century, with global warming, the glaciers around the Qinghai-Tibet Plateau have retreated completely [3], and the snow cover days and snow depth have obviously decreased [4]. Glaciers have retreated and snow has melted, exposing rocks and forming a new desert. Rock deserts affect regional differences in global climate change response through the interaction of landing surface and air. Their expansion increases the probability of local natural disasters such as soil erosion. This in turn has a significant impact on the climate, hydrology and ecological environment of the plateau, and also brings great challenges for the protection of the ecological environment and the development and utilization of land resources.

Some scholars tried to monitor desertification on the Qinghai-Tibet Plateau by remote sensing in the early 21st century. From existing monitoring research, based on the land use perspective, the desertification of the plateau was monitored by dividing the types of desertification [5,6]; based on the evaluation of ecological quality, vegetation degradation and vegetation precipitation utilization rate, the desertification of the plateau was monitored through the degree of desertification [7,8,9]. On the other hand, from the perspective of the land use and ecological quality assessment, desertification monitoring on the plateau was mainly based on the Landsat series data with medium resolution, and the monitoring results showed that the overall desertification continued to worsen with ecological degradation [10,11]. From the perspective of the vegetation degradation and vegetation precipitation utilization rate, MODIS data and NDVI products were widely used. The monitoring results show that desertification in the Qinghai Tibet Plateau presents a weak improvement trend [12]. Therefore, different remote sensing data sources would cause uncertainty in the monitoring results, using a different remote sensing monitoring technology and index system. Although some achievements have been made in monitoring desertification in the Qinghai Tibet Plateau, the attention paid to the rock desert caused by ice–snow melting is scarce. Against a background of unobvious human disturbance, dynamic changes in rock desert were a supplement to the remote sensing monitoring of ice and snow changes in the Qinghai-Tibet Plateau.

The monitoring of rock desert in ice–snow melting areas needs to be based on studies of both the ice–snow melting area and rock desert. The Landsat series can be used as a data source with high spatial resolution, radiation quantification level and time resolution. Landsat data have been providing continuous global land cover information since 1972 [13,14,15]. This is one of the main remote sensing data sources for regional snow and ice monitoring on the plateau and is also suitable for monitoring the different types of rock desert in ice–snow melting areas. Ice–snow melting areas can be monitored using snow and ice extraction methods. The S3 snow index model has high accuracy [16] and is suitable for situations when extraction accuracy is high and vegetation coverage is high [17]. Remote sensing methods to extract land desertification information can be used as a basis for rock desert extraction. Remote sensing monitoring data of the dynamic changes of ice and snow on the Qinghai-Tibet Plateau are diverse and the research methods are mature. Remote sensing monitoring of the dynamic changes of rock desert in ice–snow melting areas has not yet been a subject of a systematic study. This is main subject of this paper.

Rock desert is formed after the melting of snow and ice in the Qinghai-Tibet Plateau. This process not only leads to local ecological degradation, but more importantly, the formation of rock desert may change the regional climate at a large scale through the interaction of landing surface and air [18]. Therefore, studying the dynamic changes in the rock desert in the Qinghai-Tibet Plateau caused by climate warming and the relationship of this process with the climate are important for a more comprehensive understanding of the ecological situation on the Qinghai-Tibet Plateau. In this study, TM data from 1990 and OLI data from 2020 were taken as information sources. The snow melting areas in the last 30 years in the study area were extracted based on changes in snow lines. Remote sensing monitoring of rock desert in ice–snow melting was carried out by various methods to obtain the proportion and spatial distribution of different rock desert types. This provides support for the prevention of natural disasters such as soil erosion, debris flow, landslide, collapse and sandstorm. It is also useful information for the protection of the ecological environment, in the context of the construction of the national park of the Qinghai-Tibet Plateau.

2. Research Area and Data Source

2.1. Overview of the Study Area

Yushu Prefecture is located in the hinterland of the Qinghai-Tibet Plateau, between 31°35′–36°30′N and 89°35′–97°55′E. It is dominated by permafrost [19], and is the birthplace of China’s Yangtze River, Yellow River and Lancang-Mekong River (Figure 1). The abundant snow resources in the region play an important role in the feedback and regulation of water and energy circulation in the Qinghai-Tibet Plateau, and even the whole high Asia region. Seasonal variations of snow cover also affect snowmelt river flow, which has a significant impact on the water resources and economic activities of the relevant river basins [20]. The Kunlun Mountains in the north, and Tanggula Mountains in the south, are natural barriers between the north and south of Yushu Prefecture, and the terrain is generally low in the southeast and high in the northwest. The highest altitude is 6824 m at Bukataban Peak, and the lowest altitude is 337 m at the water surface level of Jinsha River flowing out of Qinghai. The average altitude is 4963 m. The terrain is wide, permafrost is widespread, and periglacial landform features are typical. It has a typical plateau alpine climate, with an average temperature of 0.4 °C and annual precipitation of 487.7 mm. Precipitation from May to September accounts for about 85% of the total precipitation in the whole year [21,22], and snowfall occurs every year from October to June of the following year. Yushu Prefecture has good regional representativeness and is an ideal area to study the rock desert in the ice–snow melting area of the Qinghai-Tibet Plateau.

2.2. Data Source and Preprocessing

The following three data sources within Yushu Prefecture were selected: (1) Landsat series remote sensing image data, including Landsat-5 TM from 1990 and Landsat-8 OLI from 2020, from the USGS website (http://glovis.usgs.gov/, accessed on 1 December 2021). Since there is no suitable image data for individual maps in 1990/2020, we use the images of the nearest years instead (

Table 1. The information of the dates of the Landsat images.

Path/Row	TM Image Date	OLI Image Date	Path/Row	TM Image Date	OLI Image Date
134/036	1990/08/30	2020/09/17	136/038	1989/06/22	2019/06/25
134/037	1990/08/30	2019/07/29	137/035	1990/09/04	2020/09/22
134/038	1992/09/04	2020/09/01	137/036	1990/09/04	2020/09/22
135/035	1990/08/21	2019/07/04	137/037	1990/07/02	2020/09/22
135/036	1992/07/09	2018/07/17	138/035	1991/09/14	2019/08/10
135/037	1988/08/15	2018/07/17	138/036	1991/09/30	2019/09/27
135/038	1988/08/15	2020/10/10	139/035	1990/08/17	2020/07/18
136/035	1990/08/28	2020/08/30	139/036	1989/07/29	2020/07/18
136/036	1990/08/28	2019/07/27	140/035	1991/09/28	2019/09/25
136/037	1991/09/16	2020/08/30

). To ensure the accuracy of snow extraction, three image selection criteria were adopted: time, cloud cover and snow cover. We selected images from late June to early October, with cloud cover less than 10%. The snow cover was interpreted manually and visually, and the image with the least snow cover was selected from all remote sensing images meeting the conditions of time and cloud cover. We ensured that snow cover obtained by subsequent extractions was perennial snow instead of seasonal snow. (2) Digital Elevation Model data (DEM) with a spatial resolution of 30 m, from geospatial data cloud (http://www.gscloud.cn/, accessed on 1 December 2021). (3) Data of soil texture type from the World Soil Database, at 1 km spatial resolution.

Image pre-processing included radiometric correction and image registration. Radiometric correction included radiometric calibration and atmospheric correction, whose purpose was to eliminate the errors caused by sensor instruments and the influence of atmosphere on radiation. The purpose of image registration was to make the points corresponding to the same position in space in two images correspond one to one. The remote sensing image data was pre-processed by ENVI 5.3. The steps are as follows: (1) Use the three parameters of spectral radiance value of each band, solar elevation angle and shooting time for Radiometric Calibration; (2) Use FLAASH model for atmospheric correction, according to the imaging time and central latitude of the image, the atmospheric model was determined as Mid-Latitude Summer, and the Aerosol Model was selected as Rural; (3) Image registration was performed using the image-to-image method, the 2020 images were selected as the base images, the 1990 images were warp image, and the ground control points were selected to complete the registration.

3. Classification System of Rock Desert in Ice–Snow Melting Area

The bare rock mountain in desert areas is also called “rock desert” [23]. Against the background of global climate change characterized by warming, the rock desert in the ice–snow melting area studied in this paper are areas in which ice and snow cover have receded due to natural factors, such as freeze–thaw erosion, gravity erosion, wind erosion and running water erosion, and human factors. The bedrock in the ice–snow melting area is exposed, and the surface layer of rock spalls and disintegrates, forming a landform of rock blocks and cuttings, which forms the desert in ice–snow melting areas. Rock desert can be divided into bare rock, bare stone, bare gravel, bare sand and bare soil, according to the particle size of surface weathering in the ice–snow melting area. The specific classification system is shown in

Table 2. Classification system of rock desert in ice–snow melting areas. The photos were taken in Zaduo County and Zhiduo County of Yushu Prefecture from 29 April to 4 May 2021.

Code	Category	Surface Landscape Features	Photos
1	Bare rock	After ice–snow melting, rock bedding can be seen locally in mountainous areas where bedrock is exposed, mostly in areas where snow and ice have just melted. It is distributed in the upper part of the melting area.
2	Bare stone	The clastic stage dominated by mechanical crushing of bare rock is the initial stage of bedrock weathering, and the physical weathering is the main stage, in which rocks and cuttings are formed by decomposed bedrock.
3	Bare gravel	This surface is composed of coarse sand and gravel, with no bedrock exposed, few fine particles and no soil development.
4	Bare sand	This surface is composed of aeolian sand, and is the product of further weathering of loose sediments or bedrock.
5	Bare soil	The surface layer is soil, not covered by vegetation, and is mainly composed of fine soil, with high clay content. It is mostly distributed at the bottom of the melting area.

.

4. Methods

A classification of rock desert in an ice–snow melting area should start from two aspects: extraction of ice–snow melting area and classification of rock desert in the ice–snow melting area. Two periods of snow and ice cover data were extracted using the S3 snow index model [16], and a relatively complete range of ice–snow melting area was determined by snow lines. The selection of a classification index of rock desert is related to whether it can objectively and truly reflect the spatial distribution of rock desert in the study area, and it is a key factor of classification accuracy of rock desert. In this paper, based on the field investigation of the ice–snow melting area in Yushu Prefecture, surface roughness, bare soil index, spectral reflectance and texture features of images were selected as the indicators of rock desert classification. The mature and operational multi-index factor compound analysis method, multi-index principal component analysis method, object-oriented classification method, and four combinations of the three methods were selected to classify the rock desert in the ice and snow melting area. The experimental design is represented in the block diagram shown in Figure 2.

4.1. Ice–Snow Melting Area Extraction

4.1.1. Snow and Ice Extraction Based on S3 Snow Index Model and NDVI Method

The S3 snow index model is a snow extraction method proposed by Japanese scientists Saito and Yamazaki [16]. This method is based on the bands used by the Normalized Difference Snow Index to which is added the vegetation spectrum. This takes into account the reflection characteristics of ice, snow and vegetation in different bands, effectively reducing the impact of vegetation cover on ice and snow pixel recognition, and at the same time increasing the distinction between snow and clouds, and further improving the mapping accuracy of snow cover under vegetation cover. The calculation formula of S3 is:

$$S3=\frac{NIR\ast \left(Red-SWIR\right)}{\left(NIR+Red\right)\ast \left(NIR+SWIR\right)}$$

(1)

In the formula, Red represents the reflectivity of the red band, corresponding to TM Band 3 and OLI Band 4; NIR represents the reflectivity of the near infrared band, corresponding to TM Band 4 and OLI Band 5; and SWIR represents the reflectivity of the short-wave infrared band, which corresponds to TM Band 5 and OLI Band 6. The calculation result of S3 is between −1 and 1.

According to Shimamura et al. [24], a change in vegetation coverage affects the selection of the S3 threshold in snow extraction. Generally speaking, the S3 threshold of snow cover is 0.18 in winter and spring, or in areas with low vegetation coverage, while it is generally between 0.05 and 0.18 in summer and autumn, or areas with large vegetation coverage. Being dependent on the physical properties of snow and land cover in different study areas, the threshold cannot be completely fixed. It has been found that only using the S3 threshold to extract snow will result in missing points [25,26]. These missing points are mostly caused by the thin snow thickness in the snow-covered edge area and the mixing of vegetation information. Therefore, in this study we used S3 combined with NDVI to extract snow cover. To find the best threshold using a combination of S3 and NDVI data, 100 samples were selected from the ice–snow central area and the ice–snow edge area. Based on the reflectivity of red Band4, near infrared Band5 and short-wave infrared Band6 used in the S3 snow index model, the S3 value and NDVI value were plotted (Figure 3).

When S3 ≥ 0.29, and when 0.08 ≤ S3 < 0.29 and NDVI ≥ −0.12, we can estimate that a pixel is an ice and snow pixel. The snow and ice cover results extracted by the above methods still have partial missing points and multiple points. In 2020, there were some missing points in individual spots with low snow and ice coverage and low thickness. Multiple situations involved tributaries of rivers and lakes. A man–machine interactive visual interpretation was subsequently used to revise the missing and multi-divided spots.

4.1.2. Snow Line Extraction

The lower limit of the annual snow distribution height in the hottest month is taken as the snow line of that year. The snow line observed by remote sensing images is the Seasonal Snow Line, also called Transient Snow Line. This changes with the seasons, it is low in winter and high in summer [27]. The snow line at the end of the melting season in a given year is close to the definition of snow line in glaciology, that is, the boundary between the snow-covered area all year round and the maximum melting area in summer. Above this limit, the surface is covered with snow all year round [28].

The snow line is determined by the boundary method of image feature difference between glacier accumulation area and melting area. For snow line pixels identified by the S3 + NDVI snow index model, the height value of the snow line can be obtained by superposition analysis with DEM. The DEM elevation value corresponding to each snow line pixel was assigned to the snow line pixel to obtain the snow line height. The snow line data obtained by this method are not single values; each snow line is a set of elevation data. For a specific area, the average height of all snow line pixels in the area can represent the snow line height in the area [29,30].

4.1.3. Determination of Ice–Snow Melting Area in Yushu Prefecture in the Last 30 Years

The intersection of two phases of snow line data in Yushu Prefecture is the ice–snow conservation area, and the difference set (1990–2020) is the ice–snow melting area. We used the method of layer erasure in spatial analysis in ArcGIS 10.7 to define the ice–snow covered area in 1990 as the baseline. The annual ice–snow covered area is the erasing element, and the ice–snow melting area in Yushu Prefecture for 30 years can be obtained.

4.2. Rock Desert Classification Index

4.2.1. Surface Roughness

Land surface roughness (LSR) refers to the deviation degree between the real surface and the ideal surface (the geoid) in the vertical direction. The larger the deviation, the rougher the surface, and vice versa [31]. It is an important parameter to study land surface processes such as surface runoff and soil erosion and their spatial differentiation [32,33]. Generally defined as the ratio of the surface area of the surface unit to its projected area on the horizontal plane [34], it is calculated using Formula (2).

LSR = S_DEM/S_level

(2)

Bare rock, bare stone, bare gravel, bare sand and bare soil in rock desert are the result of gradual differentiation, with gradually decreasing roughness.

4.2.2. Modified Bare Soil Index (MBI)

The MBI is a bare soil extraction method proposed by Nguyen et al. (2021) for OLI images [35]. This method mainly considers the reflection of bare soil in the near infrared (0.85–0.88 um) and shortwave infrared. Shortwave infrared 1 (SWIR1: 1.57–1.65 um) mainly reflects energy in the two bands, and shortwave infrared 2 (SWIR2: 2.11–2.29 um) mainly absorbs energy. At the same time, water absorbs most of the energy from the visible to the infrared spectrum, especially in SWIR1 and SWIR2 channels. Vegetation absorbs these SWIR wavelengths, and most of the NIR energy returns to airborne sensors after being reflected by the vegetation surface. Based on this, an MBI Formula (3):

MBI = (SWIR1 − SWIR1 − NIR)/(SWIR1 + SWIR1 + NIR) + f

(3)

In the formula, f is an additional factor with a value of 0.5, which is used to redistribute index values to negative and positive values, in which the higher positive value indicates bare soil, while the negative value indicates water body and vegetation. The factor f is only convenient to determine the range of values, so that the MBI can be used to classify vacant land threshold and other types, without affecting the ability to distinguish MBI objects. MBI values range from 0.5 to 1.5; the exposed area value is ≥0.

4.2.3. Spectral Reflectivity

Land surface reflectance is the ratio of the reflected energy of the ground object to the total incident energy. This represents the reflection ability of the object in the electromagnetic spectrum. The reflectivity of ground objects is different in different wavebands, and the properties of the ground objects can be estimated by differences in their reflectance. To select more effective bands to extract various types of rock desert in the ice–snow melting area, 100 typical sampling points were selected in the study area by constructing fishing nets to draw the spectral characteristic curves of the optimal bands (Figure 4). This shows that the reflectivity of various types of rock desert has good separability and strong regularity in the first four bands (B1–B4), and B5–B7 participated in the MBI construction, so the reflectivity values of B1–B4 were finally reserved for the extraction.

4.2.4. Textural Features

Texture refers to the spatial change of image hue as a grade function. Texture feature is an overall feature, which describes the surface properties of the scene corresponding to an image. Different from color feature, texture feature calculates multiple pixels in the region. In pattern matching, texture feature can be successfully matched, even with the conditions of local deviation. As a statistical feature, texture features often have rotation invariance and a strong resistance to noise [36]. The second-order probability statistics method can be used to obtain texture features. Second-order probability statistics uses a gray tone spatial correlation matrix to calculate texture values, which is a relative frequency matrix (the frequency of pixel values appearing in two adjacent processing windows separated by a specific distance and direction). The matrix shows the occurrence number of the relationship between a pixel and its specific neighborhood [37].

Taking the B1-B4 bands of the study area as the data sources, the co-occurrence measures tool in the texture module of ENVI 5.3 software was used to select variance, contrast, dissimilarity, information entropy, second moment and correlation. Six texture filters were used for texture analysis, and the texture processing window was 3 × 3 in steps of 1, a total of 24 texture factors were obtained.

4.3. Multi-Index Fusion Rock Desert Classification

4.3.1. Multi-Index Factor Compound Analysis (MIFCA)

Three indexes, namely surface roughness, spectral reflectance and bare soil index, were selected for compound analysis. The weight of each influencing factor was determined by the entropy method [38], the hierarchical structure diagram of influencing factors was constructed according to the relationships among the factors, and the grades were evaluated according to their importance degree, with test coefficient CR < 0.1 and passing the consistency test. The indicators and weights used in the estimation model of rock desert spatial distribution are shown in

Table 3. Index and weight of Rock desert size (RDS) calculation model.

Level-1 Indicators	Level-2 Indicators	Weight
Surface roughness	Surface roughness	0.1037
Spectral reflectance	Band1	0.1981
	Band2	0.1966
	Band3	0.1954
	Band4	0.1978
MBI	MBI	0.1082

.

The weighted summation of spatial superposition analysis was carried out by integrating the surface roughness, spectral reflectance and bare soil index. We proposed a model for estimating the spatial distribution of the rock desert, as shown in Formula (4).

RDS = ∑ W_i (LSR + Rn + MBI)

(4)

In the formula, RDS is the result of the spatial distribution of the rock desert, W_i is the sensitivity weight of the i-th factor, LSR and MBI represent the surface roughness and bare soil index, respectively, and Rn is the reflectance value of the nth band in the spectrum, with n = 1,2,…,4. Using this formula, we obtained the spatial distribution of rock deserts in Yushu Prefecture. We determined the range of various types of rock deserts based on the statistical results of typical sample data. Samples were divided into bare rock, bare stone, bare gravel, bare sand and bare soil.

4.3.2. Multi-Index Principal Component Analysis (MIPCA)

With the principal component analysis method, we selected the principal components of remote sensing and geoscience factors, and then established a regression model between the principal components and the particle size of surface material in the melting area. Because the principal components are orthogonal to each other, multiple principal components are independent from each other. This not only keeps the main characteristics of remote sensing and geoscience information, but also avoids multiple correlations among variables [39]. In this study, we therefore selected the bare soil index, spectral reflectance (B1–B4) and texture features of four image bands. Twenty-nine remote sensing and geo-science factors were used as the correlation factors of rock desert surface particle size.

The calculated sample variance of the first five principal components (PC1–PC5) reached 98.4%, reflecting the main information of the sample. Reducing the number of variables to five not only keeps the main information of the sample, but also reduces the dimension and simplifies the model. The first principal component and the second principal component had a cumulative contribution rate of 90.8%. The surface roughness was extracted from 500 typical samples selected from the rock desert. The samples were divided into 375 modeling samples and 125 verification samples. Using ArcGIS 10.7 software, the 375 sample plots in the modeling sample were added to the obtained five principal component layers, and the five principal component values of 375 sample plots were calculated. Then, through SPSS software, with the principal component value as the independent variable and surface roughness as the dependent variable, a multiple linear regression analysis was carried out. The resulting estimation model of surface gravel size in the ice–snow melting area is shown in Formula (5).

RDS = 1.171 + 0.170PC1 − 0.030PC2 − 0.373PC3 + 0.805PC4 − 0.365PC5

(5)

The R² of the regression equation was 0.694 (p < 0.01), and the test results of each coefficient were shown in

Table 4. Test of regression model coefficient.

Coefficient of Regression Equation	Coefficient Value	Standard Error	p-Values
Constant	1.171	0.010	0.000
PC1	−0.017	0.010	0.003
PC2	−0.030	0.035	0.037
PC3	−0.373	0.102	0.000
PC4	0.805	0.369	0.030
PC5	−0.365	0.342	0.026

. The constant term, PC1 coefficient and PC3 coefficient of the regression equation were all tested at the level of p < 0.01. PC2 coefficient, PC4 coefficient and PC5 coefficient passed the significance test at the level of p < 0.05. 125 samples that did not participate in the establishment of regression model were used as verification points to test the stability and prediction ability of the model; the R² of the verification sample was 0.755. Therefore, the model based on multi-index principal component analysis can be used to estimate rock desert area coverage.

4.3.3. Object-Oriented Classification (OOC)

Object-oriented classification is a processing method that integrates the structure, spectral features and geometric shape information of images [40]. This method takes the image objects generated after segmentation as the research objects and analyzes them by using the essential features of the image objects and topological relations between neighboring objects. This method enhances the spatial analysis ability of computer interpretation.

The OLI images in the ice–snow melting area only contain the information of rock desert. Due to the elimination of vegetation, ice and snow and other surface features, the spectral heterogeneity of the image has been enhanced, and the information contained in its pixels are more abundant, such as the structural characteristics, the texture information and the interconnection with neighboring features. Multi-scale image segmentation is used to set different segmentation scales for different objects in an image. Firstly, the spectral, shape and texture features of the object extracted from the image were used to establish the corresponding segmentation criteria; secondly, according to the criteria, the adjacent pixels with the same or similar spectral information were merged [41]. The segmentation effect determined the accuracy of the classification effect to a certain extent. In the specific implementation, the range of the segmentation threshold and the merging threshold was 0–100. The larger segmentation threshold, the smaller number of generated patches, and the larger segmentation threshold correspond to the poorer details of patches, and vice versa. Meanwhile, if the merging threshold was too low, some features would be misclassified, and a feature might be divided into many parts. Based on the above technical process, the segmentation threshold was set to 25 and the merging threshold was set at 80.

Using a comprehensive optimal segmentation scale, each type of sample database was established according to the spectral information, geometric and texture characteristics of rock desert, and the sample data set for extracting characteristic parameters was established according to the spatial and attributes information of characteristic sample points. Because all types of rock desert are scattered in space and have poor continuity, the extracted samples were selected based on the data of soil texture types combined with visual interpretation. At the same time, we ensured that a sufficient number of samples of various types was used. The number of samples of all types was about 1000. Finally, Object-oriented SVM was used for classification. The basic idea of classification is that the segmented object is the basic unit, combined with spectral and spatial information, by using the SVM classifier to conduct the classification process. The SVM classifier is constructed by the gaussian radial basis kernel function, where the penalty coefficient C and kernel width R are two mandatory parameters. In this paper, the C was taken as 100 and R was the reciprocal of the number of bands [42]. The bands we use for classification were B1-B4 bands, so the kernel width R was set to ¼ and the penalty parameter was set to 100 to perform the classification.

5. Processes and Results

5.1. Snow and Ice, Snow Line and Ice–Snow Melting Area

According to the results of the S3 snow index model + NDVI, the ice–snow covered areas in Yushu Prefecture in 1990 and 2020 were 2718.41 km² and 1379.42 km², respectively (Figure 5). They are mainly distributed in Bukadaban Peak in the middle of Kunlun Mountain, the birthplace of the Yellow River in the Yogurt Zonglie Basin in the northern foot of the Bayan Har Mountains, the Yangtze River source area at the junction of Zaduo County and Zhiduo County, and the Lancang River source of Gongzemuza Snow Mountain in the northern foot of the Tanggula Mountain. In 1990, the lowest value of the snow line in Yushu prefecture was 5067.63 m and the highest value of the snow line was 5636.83 m, and the average snow line value was 5283.41 m. In 2020, the lowest value was 5129.08 m, and the highest value was 5656.85 m, and the average snow line value was 5401.11 m. In the last 30 years, the average snow line in Yushu Prefecture has been above 5000 m, and the average snow line value in 2020a has risen by 117.70 m, compared with that in 1990a. According to statistics, the area of ice–snow melting in Yushu Prefecture in the past 30 years was 1451.04 km². In 2020, the snow-covered areas had different degrees of melting compared to 1990. The Yogu Zonglie Basin at the northern foot of the Bayan Har Mountains has almost completely melted. The area near the Bukadaban Peak located in the middle section of Kunlun Mountain saw a small amount of melting. The ice–snow melting area is shown in Figure 6.

5.2. Classification Results of Rock Desert in Ice–Snow Melting Area

Currently, there is no recognized method to classify rock desert types in ice–snow melting areas. In our study, we used combinations of three methods to try to obtain a high-precision classification of rock desert. This includes the fusion of multi-index factor compound analysis and multi-index principal component analysis, the fusion of multi-index factor compound analysis and object-oriented extraction, the fusion of multi-index principal component analysis and object-oriented extraction, and the fusion of three methods. In the fusion of the three methods, if the results extracted by the three methods are the same classification, the classification is determined. If the results extracted by two of the three methods are consistent, the results extracted by the two methods will be regarded as the results of this type. If the results extracted by the three methods are different, the classification will be determined by visual interpretation. In the fusion of two methods, if the results extracted by the two methods are of the same class, the class is determined. If the extraction results of the two methods are different, the classification is determined by visual interpretation.

Figure 7 shows the classification results of rock desert in ice–snow melting area based on seven methods: multi-index factor compound analysis, multi-index principal component analysis, object-oriented extraction and fusion scheme of the three methods. The total area of rock desert extracted by each method is slightly smaller than the melting area, due to patchiness in the melting area. In the process of extraction, it is easy to miss patches in the melting area, which makes the extraction area smaller than the melting area. The missing area of object-oriented classification is the largest, at 1.57 km², while the missing areas of other methods are all less than 1 km². In the extraction results of each method, the area of each type of rock desert is generally bare rock > bare stone > bare soil > bare sand > bare gravel. The area of bare rock and bare stone is about 900 km², accounting for more than 60% of the total melting area, which is higher than the sum of the other three types. The bare soil area accounts for about 20% of the total melting area.

5.3. Results of Accuracy Verification

The accuracy was evaluated by the confusion matrix, which was a common expression in the field of pattern recognition. It described the relationship between the attributes of sample and the type of classification results and is a common method to evaluate the performance of classifiers [43]. According to the confusion matrix, the accuracy evaluation indexes such as Overall Accuracy, Producer’s Accuracy and User’s Accuracy can be obtained, and Kappa Coefficient, which is used to measure the overall consistency between the reference data and the data to be evaluated. Overall Accuracy is the total number of correctly classified samples divided by the total number of samples; Producer’s Accuracy is the number of correctly classified pixels in a class divided by the number of all pixels in the reference data; User’s Accuracy is the number of correctly classified pixels in a certain category divided by the number of all pixels classified into that class; Kappa Coefficient is to establish an error matrix based on the ground survey data and image classification results, in which the matrix column represents the surface feature type of ground investigation, and the matrix row represents the classification type of image data. Producer’s Accuracy and User’s Accuracy can reflect the accuracy of different categories, and the Overall Accuracy and Kappa Coefficient can reflect the overall accuracy of the data set. This combination of indicators can comprehensively reflect the accuracy of the data set [44].

In order to verify the accuracy of the classification, we used the DJI UAV M210 and GPS field sampling in Yushu prefecture from 28 April to 6 May 2021. Sampling points were laid along both sides of the road in Yushu Prefecture, depending on traffic accessibility. The latitude, longitude and altitude of each sample point were recorded by GPS, and the land type was also recorded. We observed the nearby landscape and took photos to facilitate the follow-up work to check our results against ground object information. DJI M210 UAV is equipped with Chansi XT visible light lens, and DJI Polit flight control software is used to plan flight belt flight. The single flight shooting range of UAV was set as 300 m × 300 m, due to the obstruction of high mountains, the flight altitude of UAV was set at 400 m, the speed was 6 m/s, and the photographing mode was at equal intervals. To ensure the effectiveness of the flight, we set the heading overlap to 70% and the lateral overlap to 60%, and a total of 9 images were received. The images were taken by UAV, and we used Pix4D mapper software to complete the segmented image mosaic, and generated orthophoto for subsequent use.

5.3.1. Accuracy Verification of Ice–snow Melting Area

In the desert area, 114 sample points were collected by GPS. As snow was still falling in Yushu Prefecture during the field verification period, the sampling points in the ice–snow area were obtained on Google earth. The time of the images used was from July to August 2021, and 114 sampling points in the ice and snow area were collected. The accuracy verification results of the ice–snow melting area by confusion matrix are shown in

Table 5. Accuracy verification of ice–snow melting area extraction results.

Type	Producer’s Accuracy (Pixels)	User’s Accuracy (Pixels)	Producer’s Accuracy (%)	User’s Accuracy (%)	Overall Accuracy (%)	Kappa Coefficient
Ice and snow area	108/114	108/116	94.74	93.10	93.86	0.88
Rock desert area	106/114	106/112	92.98	94.64

.

It can be seen from Table 5 that the ice–snow melting area extracted based on S3 snow index model + NDVI method had better accuracy. Kappa coefficient reached 0.88, and producer’s accuracy and user’s accuracy were both greater than 90%.

5.3.2. Accuracy Verification of Rock Desert

Each rock desert image shot by UAVs was divided into 9 parts through the fishing net, and the type of the center point of each grid was extracted, a total of 81 samples. By GPS in the rock desert area, 114 samples were collected, a total of 195 samples per point. The accuracy verification results of the five types of rock desert through the confusion matrix are shown in

Table 6. Accuracy verification of rock desert type extraction results by seven methods.

Methods	Types of Rock Desert	Producer’s Accuracy (Pixels)	User’s Accuracy (Pixels)	Producer’s Accuracy (%)	User’s Accuracy (%)	Overall Accuracy (%)	Kappa Coefficient
MIFCA	Bare rock	32/42	32/42	76.19	76.19	76.92	0.71
	Bare stone	37/51	37/48	72.55	77.08
	Bare gravel	23/29	23/31	79.31	74.19
	Bare sand	23/31	23/32	74.19	71.88
	Bare soil	35/42	35/45	83.33	83.33
MIPCA	Bare rock	27/42	27/44	64.29	61.36	64.62	0.55
	Bare stone	29/51	29/49	56.86	59.18
	Bare gravel	20/29	20/27	68.96	74.07
	Bare sand	22/31	22/39	70.97	56.41
	Bare soil	28/42	28/36	66.67	77.78
OOC	Bare rock	34/42	34/42	80.95	80.85	79.49	0.74
	Bare stone	40/51	40/50	78.43	80.00
	Bare gravel	25/29	25/31	86.21	80.65
	Bare sand	23/31	23/31	74.19	74.19
	Bare soil	33/42	33/41	78.57	80.49
MIFCA and MIPCA	Bare rock	30/42	30/42	71.43	71.43	70.26	0.63
	Bare stone	32/51	32/45	62.75	71.11
	Bare gravel	22/29	22/34	75.85	64.71
	Bare sand	22/31	22/36	70.97	61.11
	Bare soil	31/42	31/38	73.81	81.58
MIFCA and OOC	Bare rock	36/42	36/41	85.71	87.81	83.59	0.79
	Bare stone	42/51	42/50	82.25	84.00
	Bare gravel	25/29	25/32	86.21	78.13
	Bare sand	26/31	26/34	83.87	76.47
	Bare soil	34/42	34/38	80.85	89.47
MIPCA and OOC	Bare rock	32/42	32/46	76.19	69.57	72.82	0.66
	Bare stone	32/51	32/46	62.75	69.57
	Bare gravel	24/29	24/32	72.76	75.00
	Bare sand	24/31	24/34	77.42	70.59
	Bare soil	30/42	30/37	71.43	81.08
MIFCA and MIPCA and OOC	Bare rock	34/42	34/42	80.95	80.95	78.97	0.73
	Bare stone	38/51	38/53	74.51	71.70
	Bare gravel	24/29	24/28	82.76	85.71
	Bare sand	25/31	25/35	80.65	71.43
	Bare soil	33/42	33/37	78.57	89.19

.

It can be seen from Table 6 that the accuracy of the object-oriented classification method was the highest in single method extraction, and Kappa coefficient reached 0.74. The classification accuracy of this method was high for all types of rock desert. The accuracy of bare gravel was the highest, and the producer’s accuracy and user’s accuracy were more than 80%. The texture characteristics of bare gravel were obvious. The standard pseudo color synthesis method uses a gray and white hue. The identification accuracy of bare gravel was the highest. The classification accuracy of the multi-index factor composite analysis method was the second highest. The method had strong recognition ability for bare soil, the producer’s accuracy and user’s accuracy were 83.33%. The classification accuracy of multi-index principal component analysis was the lowest, with Kappa coefficient only 0.55. This method has poor recognition ability for bare stone, and the producer’s accuracy and user’s accuracy were both less than 60%. The number of pixels wrongly divided into bare rock and bare gravel was higher.

In the multi-method fusion classification, the classification accuracy of multi-index factor composite analysis and object-oriented classification was the highest, with Kappa coefficient of 0.79 and producer’s accuracy of more than 80%, the highest among the seven methods. The overall accuracy was improved by 4.10% and Kappa coefficient increased by 0.05. The fusion of these two methods can not only utilize the spectral information and texture features of the image, but also integrate the information of each index. The overall accuracy and Kappa coefficient of the multi-index principal component analysis method were lower than the other two methods by more than 10%, and the accuracy was low. In the multi-method fusion with multi-index principal component analysis, the extraction accuracy was improved to a certain extent, but remained lower than the accuracy of the fusion of methods. The accuracy improved with the fusion of the two methods, but the precision of bare stone was still not high. In other words, in the fusion of two single methods with higher accuracy, the total accuracy can be further improved. If a single method with lower accuracy is involved in the multi-method fusion, the accuracy of the method with lower accuracy can be improved.

5.4. Composition Structure and Vertical Distribution of Rock Desert Types

According to the accuracy evaluation results, the classification accuracy of the fusion of multi-index factor composite analysis and object-oriented classification was the highest. The extraction results of different types of rock desert using this extraction method are shown in

Table 7. Area distribution of rock desert types.

	Bare Rock	Bare Stone	Bare Gravel	Bare Sand	Bare Soil	Total
Area/km2	453.62	471.48	86.43	186.68	252.42	1450.63
Proportion/%	31.27	32.50	5.96	12.87	17.40	100.00

.

The total area of the rock desert extracted by this method was 0.41 km² smaller than the ice–snow melting area. The various types of rock deserts are staggered in the ice–snow melting area, so there is no obvious characteristics in the spatial distribution, but a certain regularity in the vertical gradient. The area of bare rock and bare stone was 925.10 km², accounting for 63.77% of the ice–snow melting area, higher than the sum of the other three types, and mainly distributed on the bare ridge after the melting of ice–snow. Bare gravel accounted for the smallest proportion of 5.96%, mainly distributed in the areas below the ridge and above the foothills. Bare sand and bare soil area accounted for 30.27%, mainly distributed in the foothills and areas above. Generally speaking, in the vertical gradient, the bare rock and bare stone are in the top layer, the bare gravel is the transition layer, and the bare sand and bare soil are in the lowest layer.

6. Discussion

6.1. Ways to Improve the Accuracy of Remote Sensing Extraction of Rock Desert

In this paper, the classification scheme of multi-method combination is adopted. The classification accuracy is generally high, with the highest accuracy up to 83.59%. To further improve the accuracy in the future, we can start from improving the accuracy of a single method. For the low accuracy method based on multi-index principal component analysis, this may be due to the fact that the surface roughness is used as a dependent variable in constructing the estimation model of rock desert in ice–snow melting area. To some extent, the surface roughness reflects differences in the different types of rock desert surfaces, but it is not suitable to directly replace the real surface gravel size. If the surface roughness is taken as an index in the principal component analysis method and the surface gravel size is taken as a dependent variable to construct the estimation model, the accuracy of the rock desert may be improved. At the same time, the area of ice–snow melting area is limited, and the spatial continuity of the rock desert is poor. In this study, we used image data with spatial resolution of 30 m. According to the results of unmanned aerial vehicle shooting and field sampling, most of the study areas are mixed pixels (Figure 8). This affects the extraction accuracy to some extent. Therefore, to improve the accuracy of rock desert extraction in the future, the method can be optimized with the following three aspects: (1) Remote sensing images of higher spatial resolution can be chosen, and the proportion of mixed pixels reduced as much as possible; (2) A mixed pixel decomposition model [45] can be constructed by extracting endmembers (various land components) and estimating abundance (the proportion of each component), thus improving the extraction accuracy of various types; (3) The rock desert extraction index can be further optimized. In the future, we can consider adding indexes such as land cover topographic distribution index and rock exposure index to further optimize the rock desert extraction model.

6.2. How to Identify Rock Desert in Ice–Snow Melting Area and Rockdesert in Dry Area by Remote Sensing?

Rock desert in the ice–snow melting area and dry area are products of the climate, but their mechanisms of formation are different. The former is a warm climate rock desert, while the latter is a dry climate rock desert. Rock desert in ice–snow melting areas is bare land where permanent snow and glaciers are distributed on high mountains and melt under the action of climate warming. Rock desert in dry areas is bare land without vegetation formed by drought (mainly due to less precipitation). There are also obvious differences in material composition between the two types. Bare rock and bare stone are the main rock deserts in the ice–snow melting area. They have no obvious characteristics in spatial distribution but have obvious regularity in vertical distribution. Rock desert in dry areas is mainly bare sand, bare gravel and bare soil, distributed horizontally and in patches, mainly in inland basins, generally on both banks of river valleys in dry areas. Rock deserts in ice–snow melting areas are the result of freeze–thaw erosion, gravity erosion, wind erosion and running water erosion, while rock desert in dry areas is mainly caused by wind erosion causing dry denudation. The two kinds of soil erosion dynamics are completely different. Therefore, it is important to distinguish by remote sensing rock desert in ice–snow melting areas and rock desert in arid areas, to be able to understand the interaction between the underlying surface of different rock desert and the atmosphere and the relationship between soil erosion and bare land.

6.3. Dynamic Change of Rock Desert Response after Melting Due to Global Climate Change

In future research, revealing the response of regional land cover to global climate change resulting from dynamic changes in rock desert after snow and ice melting is an important goal of remote sensing monitoring in the snow-ice melting area of the Qinghai-Tibet Plateau [46]. In the research process, it is necessary to focus on the pattern and process of dynamic change in rock desert, analyzed by monitoring the melting of snow and ice and the dynamic change in rock desert in the Qinghai-Tibet Plateau on a long time scale. Remote sensing monitoring of dynamic changes in rock desert can reveal the mechanisms of climate effect in this area. Further studies should combine natural and socio-economic data such as meteorology, agriculture and animal husbandry, and analyze the natural and man-made factors influencing the dynamic changes of rock deserts in the ice–snow melting area of the Qinghai-Tibet Plateau. Research should further identify the respective contributions of natural and man-made factors to changes in the surface cover of the Qinghai-Tibet Plateau. Against the background of global change, such studies can contribute to building a low-carbon society and a pioneering area with high carbon absorption economy, and to provide a scientific reference for promoting regional sustainable development.

7. Conclusions

In this paper, Yushu Prefecture of Qinghai-Tibet Plateau is taken as the study area, TM data in 1990 and OLI data in 2020 are taken as information sources. Ice–snow melting area in the last 30 years are extracted based on snow lines, and various methods are used to monitor the rock desert in the ice–snow melting area by remote sensing. Based on the above analysis, we draw the following conclusions:

• According to the accuracy evaluation of various methods for extracting rock desert in the ice–snow melting area, a fast and efficient method of multi-index fusion for the classification of rock desert in the ice–snow melting area of Qinghai-Tibet Plateau is proposed. The rock desert in ice–snow melting area is extracted by seven schemes, including multi-index factor compound analysis, multi-index principal component analysis, object-oriented classification and four combinations of the three methods. The multi-index factor compound analysis and object-oriented classification compound method have the highest overall accuracy of 83.59% and Kappa coefficient of 0.79. The fusion of the two methods can not only utilize the spectral information and texture features of remote sensing images, but also integrate the information of each index. The fusion method can provide a fast and efficient reference scheme for extracting rock desert in the ice–snow melting area of Qinghai-Tibet Plateau.
• The fusion of two single methods with higher accuracy can further improve the total accuracy. If a single method with lower accuracy is involved in the multi-method fusion, the accuracy of the method with lower accuracy can be improved. In the single method classification, the object-oriented classification method has the highest accuracy, followed by the multi-index factor compound analysis method, and the multi-index principal component analysis method has the lowest accuracy. Kappa coefficients are 0.74, 0.71 and 0.55, respectively, and each method can reliably extract rock desert types. In the multi-method fusion, the multi-index factor compound analysis and the object-oriented classification compound method have the highest classification accuracy among the seven methods. Compared with the object-oriented extraction, the single method with highest accuracy, the overall accuracy is improved by 4.10%, and the Kappa coefficient is improved by 0.05. In the multi-method fusion involving principal component analysis, the classification accuracy is improved to some extent, but it is lower than that of the method with higher accuracy, especially in the fusion of the two methods, and the accuracy of bare stone is still not high.
• Yushu Prefecture of Qinghai-Tibet Plateau has a large amount of ice–snow melting, and the rock desert is widely distributed in this area. In the past 30 years, the average snow line value has risen by 117.70 m, and the ice and snow has melted by 1451.04 km², accounting for 53.78% of the ice–snow area. The ice–snow melting area accounted for 53.78% (1451.04 km²) of the ice–snow areas. The Yogu Zonglie Basin in the northern foot of Bayan Har has almost completely melted, while the area near Bukadaban Peak in the middle of Kunlun Mountain has incompletely melted. The area of bare rock and bare stone in the desert is 925.10 km², accounting for 63.77% of the total area of the ice–snow melting area. Bare sand and bare soil area account for 30.27% of the total area of the ice–snow melting area. There are different types of rock deserts in the melting area, which have no obvious characteristics in spatial distribution, but a certain regularity in vertical gradient. Bare rock and bare stone are located at the highest level, bare gravel is the transition layer, and bare sand and bare soil are located at the lowest level.