Classification of Alpine Grasslands in Cold and High Altitudes Based on Multispectral Landsat-8 Images: A Case Study in Sanjiangyuan National Park, China

Abstract

Land-use–cover change (LUCC)/vegetation cover plays a critical role in Earth system science and is a reflection of human activities and environmental changes. LUCC will affect the structure and function of ecosystems and a series of other terrestrial surface processes, such as energy exchange, water circulation, biogeochemical circulation, and vegetation productivity. Therefore, accurate LUCC mapping and vegetation cover monitoring are the bases for simulating the global carbon and hydrological cycles, studying the interactions of the land surface and climate, and assessing land degradation. Based on field GPS surveys and UAV data, with cloud-free and snow/glacier algorithms and the SVM classifier to train and model alpine grassland, the alpine grassland and LUCC were extracted by using Landsat-8 OLI satellite images in Sanjiangyuan National Park in this paper. The latest datasets of vegetation types with 30 m × 30 m spatial resolution in the three parks were prepared and formed. The classification results show that the SVM classifier could better distinguish the major land-use types, and the overall classification accuracy was very high. However, in the alpine grassland subcategories, the classification accuracies of the four typical grasslands were relatively low, especially between desert steppes and alpine meadows, and desert steppes and alpine steppes. It manifests the limitations of Landsat-8 multispectral remote sensing imageries in finer-resolution grassland classifications of high-altitude alpine mountains. The method can be utilized for other multispectral satellite imageries with the same band matching, such as Landsat 7, Landsat 9, Sentinel-2, etc. The method described in this paper can rapidly and efficiently process annual alpine grassland maps of the source areas of the Yellow River, the Yangtze River, and the Lancang River. It can provide timely and high-spatial-resolution datasets for supporting scientific decisions for the sustainable management of Sanjiangyuan National Park.

1. Introduction

Vegetation is a crucial indicator reflecting ecological changes. It has become a research hotspot, especially in sustainable development and the management of ecological systems, in recent decades [1,2]. Grasslands cover about one-third of the global land surface and are the most cultivated biome on Earth [1,3]. Currently, grassland degradation is considered as a major threat to the maintenance of ecological services, food security, and sustainable development, and jeopardizes the global effects of meeting goals and targets, such as the UN Decade on Ecosystem Restoration and Sustainable Development Goals [1,4]. Even in a world where climate change is soon halted, the global temperature rise will likely reach between 1.5 °C and 2 °C above preindustrial levels. This means that vegetation will likely face climate change effects that are substantially worse than those already experienced [3,5,6]. Compared with the accelerated global warming, human activities, e.g., intensified livestock grazing, construction of infrastructure, and unprecedented urbanization, have played more and more critical roles in changing the patterns and dynamics of vegetation on the Qinghai–Tibetan Plateau [7,8,9].

Vegetation classification is the initial basis for monitoring global vegetation dynamics and distribution patterns [10,11]. Traditional vegetation classification relies heavily on purely or semi-manual work, which needs a long period of time and intensive human, material, and financial resources [11,12,13,14]. In light of the literature, supervised remote sensing vegetation classification heavily relies on auxiliary experienced knowledge, such as field investigation, aviation aircraft images, and DEM slope and map algorithms, to fill the gaps caused by the low spatial resolution of remote sensing images [11,13,15,16], with manual delineation and semi-manual interpretation being important [16,17]. Geographic information was interpreted directly based on coarse remote sensing images by experts with field investigation knowledge, e.g., species types, vegetation structure, temporal phenology, and spatial distribution patterns [16]. However, the interpretation speed and efficiency heavily depended on the field experience of the interpreter, which was greatly influenced by the decisions of the interpreters. As there is no unified national/local vegetation plot database, vegetation plot-recording protocols, and agreed thematic outputs for classification, classified vegetation types are always confused with each other, especially in poorly sampled and plot-free regions [16,18,19]. There has been much uncertainty in the sub-class level integration and fusion of data in different formats, which are intensive, time-consuming, and inefficient [20,21]. With the progress of improved satellite sensors and band algorithms, remote sensing techniques have developed significantly in recent decades. Hyperspectral- and multispectral-resolution satellites have been applied to vegetation classification, with broad coverage, high timeliness, and the ability to make uninterrupted, long-term observations of vegetation, which play a dominant role in monitoring the dynamics of the global ecological system [11,22,23,24,25].

The majority of past studies observing grassland management from space used coarse-spatial-resolution sensors, such as the NASA Advanced Very-High-Resolution Radiometer (AVHRR) or Moderate-resolution Imaging Spectroradiometer (MODIS), or coarse-temporal-resolution sensors, such as Landsat [11,17,26,27]. Many land-cover maps at global and regional scales have been produced in recent years using remote sensing data, and the popular products include the International Geosphere Biosphere Programme (IGBP) global land-cover dataset [28], European Commission Joint Research Centre Global land cover for the year 2000 [29], University of Maryland land-cover map [30], the MODIS global land-cover products [31], and the finer-resolution global land cover [32]. However, most of the land-cover products have a coarser spatial resolution. Coarser-spatial-resolution (e.g., lower than several hundred meters) remote sensing data are not enough for catching the detailed grassland types and coverage change information [11,24,33]. Recent advances in medium-resolution data acquisition and accessibility have made Landsat-like spatial resolution remote sensing data a suitable choice for deriving finer-resolution grassland-cover maps [17,26,34,35].

Currently, the most popular algorithms for remotely sensing vegetation classification are multispectral remote sensing compositions, which capture the characteristics of land cover and acquire the classification results more easily compared with radar images and hyperspectral images. In addition, they are low in cost, have high time efficiency, and are suitable for multi-objective and multi-demand-driven vegetation classifications [12,16,18,26]. Although multispectral remote sensing images always have a broad field of view and cover a large land surface area, their temporal resolution is coarse [17,23]. Frequent cloudy weather will severely limit the production of high-quality remote sensing images by optical sensors, especially in the Sanjiangyuan region of the Qinghai–Tibetan Plateau [35,36]. As these regions have high elevation and alpine, rainy climates, most of the remote sensing images are more or less affected by clouds and snow, which seriously affect the identification and classification algorithms of subsequent images [11,36,37]. In addition, these regions have high elevations and are commonly frigid during the year. Accordingly, the vegetation growing season is extremely short. Cloud-free and snow-free remote sensing images are relatively scarce in the growing season [35]. Snow coverage and glaciers are the other crucial influencing factors for vegetation classification accompanied by frequent clouds in these regions [38,39]. How to effectively remove the influence of cloud and snow on images is particularly critical in image processing and application. Therefore, minimizing or removing the influence of cloud and snow is the basis for the accurate interpretation of remote sensing images and is a meaningful way to make remote sensing data reasonable and effective [38,39,40].

Cloud-removal algorithms for multispectral remote sensing images are the challenges focused upon in alpine regions at present [17,39]. Many thin-cloud-removal algorithms based on different principles and techniques have been developed and have different performances, and have been confined to modular environments in recent years [20,32,35,41]. The representative method is the Automated Cloud Cover Assessment (ACCA) system [42]. It uses multi-spectral filters and the thermal infrared band for detecting cloud and shadow, and has historically been used for cloud filtering of Landsat images. Another typical algorithm is the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) for generating an internal cloud mask [43]. For low- and mid-latitudes, LEDAPS has good performance in atmospheric correction and cloud detection. At present, the operational algorithm C-function of mask (CFmask) is used for generating feature masks in the Landsat images provided by the U.S. Geological Survey (USGS, https://earthexplorer.usgs.gov/ accessed on 20 July 2022) and estimating the pixel fraction of cloud cover in a single scene accompanied by the Landsat level 1/2 imageries [44]. The Fmask algorithm was originally designed to generate cloud, shadow, and snow masks for Landsat 4–7, and was later expanded to Landsat-8 as CFmask by utilizing the cirrus band [45]. The CFmask algorithm has been demonstrated to perform well in high altitudes, but overestimates cloud cover over snow/ice and water biomes and bright targets, and underestimates cloud cover over desert/hot biomes, primarily because it uses a fixed-scene-based threshold for all of the pixels [46]. Zhao et al. explored a cloud-removal algorithm using a low-pass filter to filter the noise and extract the low-frequency components, selected thin clouds in the Three Gorges area of the Yangtze River in eastern Sichuan, and obtained empirical results that the method is particularly suitable for removing thin clouds in large areas [47]. Li et al. proposed a novel algorithm for cloud and shadow detection based on multispectral images using a pure-image-analysis method, and applied the detection algorithm to effectively remove clouds and shadows in a data fusion framework [48]. In 2005, Wang et al. explored a high-pass wave de-clouding method based on multispectral images, which was easy to realize and could directly process cloudy images. Furthermore, this method obtained promising results for cloud-covered ETM images, and also had good performance, particularly for removing thin clouds from a large area [37]. Liu et al. analyzed the removal of thin clouds on multispectral images and proposed to select suitable de-clouding products for different atmospheric conditions and sensors; however, the heavy snow was difficult to treat [49]. Ma et al. proposed a new multispectral-based MODIS image de-clouding algorithm, but the requirements for alignment between images and data were stringent [50]. However, traditional cloud-removal algorithms do not easily remove clouds from high-elevation mountainous area optical remote sensing images. Shen et al. introduced an air-domain-filtering method based on the traditional homomorphic filtering, which has largely improved the ability to remove the clouds from images of mountainous areas [51]. In 2014, Long et al. applied the DCP algorithm for cloud removal, and used a low-pass Gaussian filter to smooth the calculated atmospheric reflection to avoid the halo artifacts generated by the DCP algorithm [52]. Subsequently, DCP has been continuously improved and used for thin cloud removal [53]. A GRS-HTM algorithm based on the improvement of the HTM algorithm was proposed [54]. It is based on the initial HTM calculation and suppresses the local area radiation centered on the ground edge points to calculate a more accurate cloud distribution map and avoid the excessive defogging phenomenon. In 2019, Xu et al. proposed a noise-adjusted principal component transformation-based cloud-removal algorithm, CR-NAPCT [55]. They proved that, the more clouds in remote sensing images, the higher the signal-to-noise ratio (SNR) of high-altitude correlation, and, conversely, the fewer clouds in remote sensing images, the lower the SNR, effectively discriminating the presence of clouds in pixels. At present, cloud-detection algorithms mainly focus on four typical categories [56]: (1) Physical rule-based algorithms based on the physical properties of clouds/cloud shadows that determine the optimal thresholds to achieve masking [42,44]; (2) Temporal change-based algorithms that compare the temporal changes in clear-sky surfaces and the sudden changes in surface reflectance caused by clouds/cloud shadows [56,57,58]; (3) Variational model-based algorithms, which construct a variational model with a priori constraints based on the prior knowledge of cloud-cover components in an image and the cloud-free image components, and detect clouds through an optimal solution of the variational model [59,60]; and (4) Machine-learning-based algorithms, involving constructing a suitable classification model, e.g., CNN, SVM, or RF, and iteratively optimizing the model parameters based on largescale training data, resulting in the model having certain generalized application capabilities [61,62].

Although a number of cloud-masking algorithms have been developed for optical sensors, very few studies have carried out a quantitative intercomparison of state-of-the-art methods in this domain. A representative work was carried out by the Cloud Masking Intercomparison eXercise (CMIX, https://calvalportal.ceos.org/cmix accessed on 20 July 2022) conducted within the Committee Earth Observation Satellites (CEOS) Working Group on Calibration and Validation (WGCV). CEOS is the forum for space agency coordination and cooperation on Earth observations, with activities organized under working groups. CMIX, as one such activity, is an international collaborative effort aimed at intercomparing cloud-detection algorithms for moderate-spatial-resolution (10–30 m) spaceborne optical sensors [46]. CMIX-I has evaluated the most representative 10 algorithms, namely ATCOR, CD-FCNN, Fmask 4.0 CCA, FORCE, IdePix, InterSSIM, LaSRC, MAJA, S2cloudless, and Sen2Cor, developed by 9 teams from 14 different organizations representing universities, research centers, and industries, as well as space agencies (CNES, ESA, DLR, and NASA). However, the CMIX-I final report [63] demonstrates that there is not an optimal cloud-detection algorithm for general utilization in the Landsat image series, as the performance of algorithms varies depending on the reference dataset [46]. The same conclusions were summarized by Tarrio et al. [64] after examining the relative performance of five different cloud-masking algorithms (Sen2Cor, MAJA, LaSRC, Fmask, and Tmask). They recommended that no one algorithm produces the best results for detecting both clouds and shadows in Sentinel-2 images. They suggested that combining the results from multiple algorithms can improve the overall accuracy of the results, but clearly require more high-quality datasets of annotated images and computational effort. Until now, there has not been a globally available cloud/shadow-masking dataset. The most influential Landsat-8 Cloud Cover Assessment Validation Data of USGU (https://landsat.usgs.gov/landsat-8-cloud-cover-assessment-validation-data accessed on 20 July 2022) are based on Fmask with C program (https://www.usgs.gov/landsat-missions/cfmask-algorithm accessed on 20 July 2022) and delivered by NASA with Landsat images. However, the fractional snow-covered area (fSCA) product is only available in Alaska and the Western U.S.—there are no data for Asia and the Tibetan Plateau (https://www.usgs.gov/landsat-missions/landsat-fractional-snow-covered-area-science-products accessed on 20 July 2022). Many researchers propose improving the already-developed detectors with large datasets of annotated samples [45,60], or developing novel, high-efficiency, and high-accuracy cloud/shadow-masking algorithms [46,56,62]. Additionally, in high-altitude mountainous areas, the number of low cloud/snow percentage images is very limited due to the short growing season—these methods have limited performance in high elevations [10,38,65,66,67]. There is a high need to develop novel, high-accuracy, and high-efficiency cloud/shadow-masking algorithms in these areas.

The Sanjiangyuan Region (SJY) is 363,000 km² (31°39′–36°12′N, 89°45′–102°23′E), which accounts for half of the area of Qinghai Province. The SJY is called the “Chinese Water Tower” due to its abundant water resources, and its name means “Three River Source”, being the source of the three most important rivers in China: the Yangtze, the Yellow, and the Lancang Rivers [68]. Sanjiangyuan National Park (SNP) includes three early established nature reserves, including the Yellow River Source Park, the Yangtze River Source Park, and the Lancang River Source Park [69]. SNP plays a critical role in maintaining the ecological balance of the Qinghai–Tibet Plateau [35]. It has a long, frigid winter, but a very short, temperate summer. As the vegetation growing season for alpine grassland is very short, there are very few satellite images for cloud/snow removal and vegetation classification, which seriously affects the recognition and classification of subsequent images. There is an urgent need for new, high-resolution data products to monitor and classify the alpine grassland in SNP. The paper takes SNP as a case study to classify alpine grasslands in cold and high altitudes based on multispectral Landsat-8 images, aiming to facilitate the methods for similar alpine regions.

2. Materials and Methods

2.1. Data Source and Pre-Processing

2.1.1. Landsat-8 Remote Sensing Images

This paper is based on the image data of the Landsat-8 Operational Land Imager (OLI) of the United States Geological Survey (USGS) (http://glovis.usgs.gov accessed on 20 July 2022) for vegetation classification and extraction. The Landsat-8 satellite can achieve global coverage every 16 days, with an imaging width of 185 km × 185 km. It maintains a good continuity with Landsat 1–7 in terms of spatial resolution and spectral characteristics. It has 11 bands, band 1 to band 7, and band 9 to band 11 having a spatial resolution of 30 m, and band 8 is a panchromatic band with a resolution of 15 m. The OLI sensor includes all the wavebands of the ETM sensor, with the major change being the adjustment of band 5 to 0.845–0.885 μm to exclude the effect of water vapor absorption at 0.825 μm. Panchromatic band 8 has a narrower range to better distinguish between vegetated and non-vegetated areas. The new blue band 1 (0.433–0.453 μm) is used for coastal zone observation, and short-wave infrared band 9 (1.360–1.390 μm) is used for cirrus cloud detection [23]. Glaciers and snow could be heavily influenced and mixed by clouds in summer rainy seasons. Therefore, a window period with the least cloud and snow cover is the optimal time. We have demonstrated that the best period is between late September and early November.

Based on the objective of vegetation classification, multi-temporal Landsat-8 remote sensing images captured in the growing seasons from July to earlier November between 2018 and 2021 in the SJY were used in this case study. The images from 2021, with a few scattered, thin clouds, thick clouds, or snow, were used as the target images. The multi-temporal images between 2018 and 2021 with fewer clouds and snow were used as the reference images for regression analysis (Figure 1).

The purpose of pre-processing multi-temporal data is to match the reference image with the target image in terms of both the spatial location and spectral characteristics, and to prepare the data for subsequent cloud area-detection and cloud removal. Landsat-8 data, similar to other TM/ETM+ data, are labeled L1T with systematic radiometric correction and geometric correction for terrain participation when released, and can be used directly without geometric correction. According to the objective of remote sensing monitoring applications, the images were radiometrically calibrated and atmospherically corrected to eliminate the influence of the sensors themselves and factors such as the scattering of atmospheric molecules and aerosols in images of different phases on the reflections of features, to obtain the true reflections of grassland features, and then use their rich wavelength spectral information.

(1) 

Radiometric Calibration. Radiometric calibration converts the DN value of the raw image to the reflectance at the top of the atmosphere, which eliminates the response differences between different sensors in the same image, and also eliminates the effect of solar altitude angle on the image. The process uses radiometric calibration in ENVI 5.1 software to process the raw image, where the associated irradiance conversion parameters are read from the image head files for FLAASH atmospheric correction;

(2) 

Atmospheric Correction. All Landsat-8 data are atmospherically corrected by the FLAASH atmospheric correction module in ENVI 5.1. FLAASH uses the MODTRAN4+ radiative transfer model to effectively remove water vapor/aerosol scattering effects. At the same time, based on pixel-level correction, the “proximity effect” of the cross-radiation between the target pixel and adjacent pixels is corrected, and high-precision ground-object-reflectivity data can be obtained.

2.1.2. GPS and UAV Sample Vegetation Types

In order to finely extract the vegetation types in the SNP, this research group conducted an intensive field investigation in the three parks in August 2018. Typical vegetation types were sampled by GPS positioning (Figure 2), and a ground-verification-point dataset for the three parks was developed (http://sjynp.tpdc.ac.cn/zh-hans/ accessed on 20 July 2022). Furthermore, a visible-light camera equipped with DJI Phantom 4 drone was utilized to photograph the vegetation types in some typical vegetation areas. Ultimately, the ground-vegetation-type information datasets of the SNP were established (http://sjynp.tpdc.ac.cn/zh-hans/ accessed on 20 July 2022). After careful examination and exact classification, 1/3–1/2 were reserved as verification samples, and the rest were treated as training samples for the classification modeling.

2.2. Cloud-Removal Algorithms

Multispectral remote sensing images are popular data sources for remote sensing retrieval and remote sensing monitoring of the environment on a global or regional scale. Due to climatic influences, frequently cloudy weather exacerbates difficulties for optical sensors to obtain high-quality remote sensing images, especially in the vast QTP. It is a challenge to obtain completely cloud-free remote sensing images. Most remote sensing images from high elevations are more or less affected by cloud or snow when they are acquired. Therefore, how to remove the influences of cloud and snow has been considered a complex issue in image processing and classification.

For practical applications, such as large-scale vegetation remote sensing retrieval and dynamic monitoring, it is not only necessary to effectively recover the image areas obscured by thin clouds, thick clouds, and cloud shadows, but also to remove or diminish the influence of clouds and shadows with the minimum loss of image information, and to maintain the qualities of images after de-clouding. This process is the prerequisite for practical utilizations, and it provides high-quality pre-processing products for remote sensing retrieval and monitoring.

(1) 

Cloud detection

From the spectral digital number (DN) values obtained from the sampling paths (Figure 3), the cloud area could be directly identified through the cirrus band (Band 9) of Landsat-8 OLI. When the empirical representative spectral value of Band 9 is greater or equal to a specific threshold Δ_Cirrus, the cloud area could be directly calculated. The DN values of other features in Band 9 were always less than this threshold. With this condition for cloud area-detection, the cloud area-detection model of Landsat-8 was obtained as Equation (1):

DN_T9 ≥ Δ_Cirrus

(1)

where DN_T9 is the spectral value of the Landsat-8 Band 9 image, and Δ_Cirrus is the cloud threshold in the cirrus band. The empirical cloud threshold in this study was 5400.

(2) 

Shadow detection

The reflectance of the shaded area and the water body showed a rapid decrease from Band 2 to Band 7 (Figure 3), and the reflection values were significantly less than those of other feature types. If the differences of each band between the objective imagery and referenced imagery are calculated, the mean absolute difference M_abs must be greater than the threshold ΔSM. In addition, the reflection values of the cloud-shaded areas in bands 2 to 7 were all close to each other, while the reflections of water bodies in bands 2, 3, and 4 were slightly higher than the values in bands 5, 6, and 7; therefore, water bodies could be removed from the areas with lower reflectance values. By using this condition for detecting cloud shadow areas, the cloud shadow area detection model of Landsat-8 OLI is (Equation (2)):

M_abs × Sign(N − n + 1) ≥ ∆SM

(2)

where n is the number of bands, M_abs is the mean absolute difference of 2 images in two periods, ∆SM is the cloud shadow area detection model, N is the cloud shadow area discriminant function, Sign() is the sign function, and the values are 1, 0, and −1 according to the positive, zero, and negative of the variables, respectively. For water bodies, which are easily mixed with shadows, the normalized water body index NDWI (Equation (3)) has been used to eliminate them from shadows [36]:

$$NDWI=\frac{T_3-T_5}{T_3+T_5}$$

(3)

where T₃ and T₅ are the Band 3 (green) and Band 5 (near-infrared) of the referenced image.

(3) 

Cloud and shadow removal

The cloud-free area in the target image and the corresponding area in the referenced image were subjected to spectral linear regression analysis to obtain the linear regression coefficients between the two images. Since the temporal phases of the two Landsat-8 images were very close, the correlation coefficients between each band were relatively high. The high correlation coefficients not only indicated the high spectral similarity between the corresponding bands of the two images, but they also demonstrate that the variations in the ground coverage are very slight. The reference image after atmospheric correction could objectively represent the cloud-cover target image. Therefore, in order to better eliminate the differences in spectral features between the target image and the reference image, the linear regression coefficients between the target image and the reference image were calculated to establish the linear fitting models by the least-squares method. Finally, Landsat-8 cloud-free products were obtained (Figure 4).

2.3. Snow- and Glacier-Removal Algorithms

Cloud contamination can significantly limit the signal quality of snow/glacier property detections made by satellite optical–infrared spectrum remote sensing. For multi-band remote sensing images, the spectral–band ratio method has proved to be a simple, highly efficient, robust, and accurate technique to extract glacier outlines [70,71,72,73,74]. The slight differences in the spectral features of the objects are enhanced, which helps to distinguish the types of objects. The method avoids the problems of sensor saturation and shadowed areas, and discriminates debris-mantled ice and ice-marginal water bodies [73]. We used Landsat spectral bands to calculate the band ratio (Equation (4)):

Ratio = CHn/CHm

(4)

where n is the band number of the red spectral (Red) band (Band 4 in OLI) or the near-infrared spectral (NIR) band (Band 5 in OLI), and m is the band number of short-wave infrared spectral (SWIR 1) band (Band 6 in OLI).

To eliminate snow-cover influences, the normalized-difference snow index (NDSI) was used to distinguish the snow zones (Equation (5)):

NDSI = (CHn − CHm)/(CHn + CHm)

(5)

where n is the band number of the visible spectral band, and m is the band number of the near-infrared spectral band, e.g., Band 3 and Band 6 in OLI. This is based on the difference between the strong reflection of visible radiation and the near-total absorption of shortwave infrared wavelengths by snow [74,75]. The NDSI has been effective in distinguishing snow from similarly bright soil, vegetation, and rock, as well as from clouds in Landsat imagery [72]. Although the spectral band ratio method could quickly extract glacier outlines, it needs several attempts to determine the threshold of the ratio. Inappropriate thresholds will cause misclassification between glaciers, bare land, and water bodies. On-screen digitizing by the manual delineation of glacier ice is time-consuming and labor-intensive; however, it could make up for the shortcomings of the spectral band ratio method. It is still widely used in combination with supervised classification techniques, especially when the analyst is knowledgeable about snow/glaciers [71,74,75,76]. The DEM is employed to eliminate the effects of perspective distortion and to reduce the topographic effect of remotely sensed data [77,78]. Combining the above-mentioned methods, the outline of the glacier and snow can be distinguished from clouds and other objects by the spectral band algebraic operation.

In this study, pseudocolor synthesis was performed on the 7, 5, and 3 bands of Landsat-8 cloud-free images produced by the de-clouding algorithm, and the training samples were selected sequentially based on the ground sampling point dataset for vegetation recognition by the support vector machines (SVM) classifier [79,80]. The SVM algorithm is one of the most accurate and robust algorithms in data mining and can be used to classify linear and nonlinear data. In most cases, SVM-based classifiers can achieve better classification accuracy than other widely used classifier techniques [81]. The SVM algorithm was first proposed by Cortes and Vapnik [82] to classify data that were linearly separable and later generalized to nonlinear states. By using this algorithm, the data were categorized into two phases of train and test sets, and, to validate the data, cross-validation techniques, such as k-fold, holdout, or leave-n-out training methods, were applied [82,83]. In order to find the optimal value of the penalty parameter (C) and the kernel parameter (σ), the holdout cross-validation was applied to perform the classification with the highest accuracy. The SVM classifier outperformed univariate decision trees, maximum likelihood, and backpropagation neural network classifiers, and with a limited number of training samples, there was almost no Hughes phenomenon [79,80,81,84]. There are four main advantages of SVM. Firstly, it has a regularization parameter, which makes the user think about how to avoid over-fitting. Secondly, it uses the kernel trick, so one can build on expert knowledge about the problem by engineering the kernel. Thirdly, SVM is defined by a convex optimization problem (no local minima) for which there are efficient methods. Lastly, it is an approximation to a bound on the test error rate, and there is a substantial theoretical knowledge base underlying the model [81]. Owing to these advantages, SVM is widely used in remote sensing classification [80,81,84,85]. In this paper, 1/3 of the samples were considered as test data each time, and the remaining 2/3 were taken as training data. The common radial basis function (RBF) is used in SVM classification modeling. The overall technical flowchart is shown in Figure 5.

2.4. Vegetation Classification Tree

This study utilizes the currently widely used IGBP DISCover land-use type classification system of the IGBP (International Geosphere–Biosphere Programme). The system adopts the USGS classification method, and it was furtherly developed by the IGBP [28,86]. According to the system, there are 17 classification types on Earth (

Table 1. IGBP DISCover land cover classification definition system.

Code	Classification	Description
1	Evergreen Needleleaf Forest	Land dominated by trees with a percent canopy cover of >60% and height exceeding 2 m. Almost all trees remain green all year. Canopy is never without green foliage.
2	Evergreen Broadleaf Forest	Land dominated by trees with a percent canopy cover of >60% and height exceeding 2 m. Almost all trees remain green all year. Canopy is never without green foliage.
3	Deciduous Needleleaf Forest	Land dominated by trees with a percent canopy cover of >60% and height exceeding 2 m. Consists of seasonal needleleaf tree communities with an annual cycle of leaf-on and leaf-off periods.
4	Deciduous Broadleaf Forest	Land dominated by trees with a percent canopy cover of ≥60% and height exceeding 2 m. Consists of seasonal broadleaf tree communities with an annual cycle of leaf-on and leaf-off periods.
5	Mixed Forest	Land dominated by trees with a percent canopy cover of >60% and height exceeding 2 m. Consists of tree communities with interspersed mixtures or mosaics of the other four forest cover types. None of the forest types exceeds 60% of the landscape.
6	Closed Shrublands	Lands with woody vegetation less than 2 m tall and with shrub canopy cover of >60%. The shrub foliage can be either evergreen or deciduous.
7	Open Shrublands	Lands with woody vegetation less than 2 m tall and with shrub canopy cover between 10–60%. The shrub foliage can be either evergreen or deciduous.
8	Woody Savannas	Lands with herbaceous and other understorey systems and with forest canopy between 30 and 60%. The forest cover height exceeds 2 m.
9	Savannas	Lands with herbaceous and other understorey systems and with forest canopy of 10–30%. The forest cover height exceeds 2 m.
10	Grasslands	Lands with herbaceous types of cover. Tree and shrub cover is less than 10%.
11	Permanent Wetlands	Lands with a permanent mixture of water and herbaceous or woody vegetation that cover extensive areas. The vegetation can be present in either salt, brackish, or fresh water.
12	Croplands	Lands covered with temporary crops followed by harvest and a bare soil period (e.g., single and multiple cropping systems). Note that perennial woody crops will be classified as the appropriate forest or shrub land-cover types.
13	Urban and Built-Up	Land covered by buildings and other man-made structures. Note that this class will not be mapped from the AVHRR imagery, but will be developed from the populated places layer that is part of the Digital Chart of the World (Danko 1992).
14	Cropland/Natural Vegetation Mosaic	Land with a mosaic of croplands, forest, shrublands, and grasslands in which no one component comprises more than 60% of the landscape.
15	Snow and Ice	Land under snow and/or ice cover throughout the year.
16	Barren or Sparsely Vegetated	Land of exposed soil, sand, rocks, or snow that never has more than 10% vegetated cover during any time of the year.
17	Water Bodies	Oceans, seas, lakes, reservoirs, and rivers. Can be either fresh or salt water.

).

Based on the unique characteristics of the QTP, such as high altitude, severe cold, and short vegetation growing season, the vegetation and land-cover types in the SJY were classified into 7 major types based on the China Vegetation Type System 2007 [87]. Firstly, based on the multi-year NDVI average, the areas with NDVI less than 0.05 were judged as unvegetated areas [88,89]. The alpine grassland types were subdivided into four categories: swamp meadow, alpine meadow, alpine steppe, and desert steppe. This division can objectively and reasonably reflect the alpine vegetation types and their characteristics in the SJY. It can also achieve a more refined classification of vegetation types within the capacity of multispectral remote sensing. Considering that SNP is a national park area with very little agricultural, settlement, or construction land, it is difficult to interpret these plots into a separate land-use type alone in the overall system. In the main categories, non-vegetated sandy land, mudflats, saline land, and other types are combined into bare land. The division scheme can objectively reflect the vegetation characteristics of the SJY and realize the reasonable classifications of subcategories in the context of ecology and hydrology. The overall vegetation classification system and remotely sensed feature characteristics are shown in

Table 2. Vegetation classification system and field characteristics of Sanjiangyuan National Park.

Code	Name	Illustration	Description
1	Water		The color composite map is black, homogeneous, and has contiguous blocks with clear boundaries, which are easy to distinguish from other grasslands.
2	Swamp meadow		The formation of richly colored patches of standing water and low-lying areas near rivers and lakes with seasonal standing water, interspersed with puddles of standing water, is extremely easy to identify in combination with the topography.
3	Alpine meadow		In the upper part of the marsh meadow, the color is dark green, clearly distinguished from other vegetation types, and interspersed with alpine steppe vegetation types.
4	Alpine steppe		The color is shiny and light green, and interspersed with alpine meadows and desert steppes, which are more difficult to distinguish.
5	Desert steppe		Near the bare ground, the color can be clearly distinguished from bare soil, with sparse vegetation and bright colors.
6	Glacier/snow		Strong spectral albedo, easy to distinguish from other features, bright blue or white distribution, and obvious morphological features combined with the topography, making it easier to identify.
7	Bare land		Lack of vegetation cover, dark purple or dark red color, and obvious textures of wind-forming and water-forming effects.

.

3. Results and Discussion

3.1. Classification Accuracy and Validation

Based on the SVM classifier, this paper classified the grassland based on Landsat-8 multitemporal cloud-free remote sensing images of the three parks in the SJY. A land cover/vegetation type dataset (2021) of the three parks in SNP by manual verification and correction of misclassification was produced. As the classification accuracy was the main issue for its application, we estimated the errors of our SVM method in this study. The average minimum distance matrix (Least to Most) between each classified type was calculated based on the statistics of the classifier sample points (Figure 6).

According to the statistical results, the SVM classifier can better distinguish the four main grassland types, and the overall accuracy of classification is relatively high. The average Kappa coefficient among the four grassland types in the three parks is 0.8366, the overall classification mapping accuracy can reach 84.52%, and the mapping accuracy for users could reach as high as 87.67%. The statistical results indicate that the average minimum distance between some first-level classes reached the maximum value of 2.0. For example, the differences between bare land/water bodies and other classes were relatively high, and the classification accuracies were high overall. However, within alpine grasslands, the classification accuracies of the four sub-grasslands were relatively low, especially between desert steppe and alpine meadow (1.3297), and desert steppe and alpine steppe (1.4171). We believed that this was highly related to the limited number of ground GPS sample points and the low heterogeneity of the selected sample points [84,90,91]. The high degree of similarity of the spectral characteristics of the four grassland types was also considered to be the main reason [23,81,85]. On the other hand, the low classification accuracy has demonstrated the limitations of Landsat-8 multispectral remote sensing imageries in finer-resolution grassland classifications of high-altitude alpine mountains [23,36,92].

In order to compare the accuracy of the produced datasets with other similar data products, the MODIS 500 m vegetation classification dataset MCD12Q1 (https://lpdaac.usgs.gov/products/mcd12q1v006/ accessed on 20 July 2022), ESA-CCI Global Land Cover 300 m (http://maps.elie.ucl.ac.be/CCI/viewer/download.php accessed on 20 July 2022), and Global Land Cover 30 m (GLB30, http://www.globallandcover.com/ accessed on 20 July 2022) were then selected for the comparison of classification accuracy in this paper.

From the two randomly selected validation sites in the Yellow River Source Park (Figure 7), the data produced in this paper have outstanding classification advantages in terms of the type richness of land use/vegetation cover and the spatial and temporal resolutions of the dataset. The performances of these datasets demonstrate that the classification results of this paper can optimally distinguish the four alpine grassland types. However, they could not be classified or performed for the other multispectral remote sensing datasets due to their low spatial resolutions or coarse classification systems. In these datasets, for example, GLB30 had the most consistent spatial resolution with the dataset in this paper. It only classified alpine grassland types as “grassland”, and classified swamp meadow as “wetland”. The classification schedule cannot reflect the local characteristics of SNP in a finer resolution. The comparisons demonstrate that the dataset of this paper has the advantages of more alpine grassland types and higher spatial and temporal resolutions.

3.2. Classification Results

(1) 

Yellow River Source Park

The Yellow River Source Park is located in the eastern part of the SNP (Figure 8), with a total area of about 19,000 km². It has the top-two largest plateau freshwater lakes, Zhaling Lake and Eling Lake, with a total water area of 1524.6 km², in the territory, accounting for 8% of the total park area at the source of the Yellow River. On the whole, the north is low, with relatively little precipitation and mainly desert steppe areas, while the southern part is relatively high, with abundant precipitation and mainly alpine meadow and alpine steppe areas. According to the classification results (Figure 8), the most widely distributed areas were alpine steppe and desert steppe. The areas were 6189.9 km² and 6405.9 km², respectively, accounting for 32.6% and 33.7% of the total area of the park. The dominant species of the alpine steppe were Elymus nutans, Roegneria nutans, Kobresia tibetica, and Carex spp. Swamp meadows were mainly distributed near lakes or the Yellow River channel on flat terrain, with an area of 1326.1 km², accounting for 7% of the park’s total area. Swamp meadows in this area were dominated by Kobresia littledalei, Kobresia tibetica, and Kobresia prattii. Alpine meadows were mainly concentrated in the high-altitude mountains in the south, with an area of 2175.2 km², accounting for 11.4% of the total area. The species of the alpine meadow were dominated by Kobresia pygmaea, K. humilis, K. setchwanensis, K. capillifolia, Polygonum sphaerostachyum, Kobresia tibetica, Carex lanceolata, C. muliensis, and C. meyeriana. In addition, mudflat-dominated bare land was distributed near the Yellow River channel, accounting for 6.5% of the park area. There were glaciers or permanent snow scattered on the tops of high mountains, but these covered less than 1% of the total area.

(2) 

Yangtze River Source Park

The Yangtze River Source Park has the largest area among the three parks of SNP, with a total area of about 90,500 km². It is located in the northwest of the SNP (Figure 9). The elevation of the eastern part of the SJY is generally low and gradually increases from east to northwest. The Yangtze River Source Park has the highest average elevation among the three parks of SNP and the same elevation distribution trend as SJY. According to the classification results (Figure 9), the whole park was dominated by alpine vegetation. It was a typical alpine meadow, alpine steppe, and desert steppe area. The north and west were mainly dominated by desert steppes and bare lands due to the high altitudes and abundant glaciers and water bodies. The desert steppes covered 19,776 km² and bare lands covered 17,230 km², accounting for 21.7% and 18.9% of the park’s area, respectively. The dominant species of the desert steppe were Carex moorcroftii, Corex pseudofoetida, Stipa purpurea, Festuca ovina, Artemisia spp., Ceratoides compacta, Ceratoides latens, Stipa glareosa, Ephedra gerardiana, Orinus thoroldii, and Pennisetum flaccidum. The central and eastern areas with lower elevations were mainly alpine steppes and swamp meadows. The areas were 26,294 km² and 14,240 km², accounting for 28.9% and 15.6% of the total area of the park, respectively. The dominant species of alpine steppe and swamp meadow of this area were Kobresia tibetica, Kobresia prattii, Kobresia littledalei, Kobresia pygmaea, Kobresia deasyi, K. humilis, Carex lanceolata, C. muliensis, and C. meyeriana. In general, the vegetation in the Yangtze River Source Park has a very short growing season, very little above-ground biomass, relatively weak grassland-carrying capacity, and high vulnerability in alpine grassland ecosystems.

(3) 

Lancang River Source Park

This area is located in the south of SJY and is the smallest of the SNP. It is influenced by warm and humid moisture from the southeast Asia monsoon. The annual mean precipitation is about 400–600 mm, which is relatively high for this region. Based on the classification results in this paper (Figure 10), the area was mainly dominated by alpine steppe, alpine meadow, and swamp meadow and was a typical plateau alpine grassland area (Figure 10). Alpine grassland covered 3207.7 km², alpine meadow covered 2877.7 km², and swamp meadow covered 2241.6 km², accounting for 23.5%, 21.1%, and 16.4% of the total area of the park, respectively, and 61% of the total area of the park. The alpine steppe was mainly dominated by Kobresia setchwanensis, Kobresia tibetica, Kobresia deasyi, Kobresia prattii, Kobresia littledalei, Kobresia pygmaea, Kobresia humilis, Carex lanceolata, and Corex spp. The alpine meadow was mainly dominated by Elymus mutans, Roegneria nutans, Corex spp., Kobresia pygmaea, Kobresia setchwanensis, Kobresia littledalei, and Thylacospermum. The elevation of the western part of the region was relatively high, with more glaciers and bare land distribution, while the northern part and river mudflats were distributed with desert steppes, dominated by sparse vegetation, such as Sanssurea, Stipa aliena, Stipa purpurea, Stipa breviflora, Stipa glareosa, Festuca ovina, S. purpurea var. arenosa, Littledalea racemose, Carex moorcroftii, Ceratoides compacta, Ceratoides latens, Salsola abrotanoides, Artemisia sphaerocephala, Puccinellia spp., Polygonum sibiricum, and Orinus thoroldii. In addition, there were significantly developed networks of braided water systems, abundant water bodies, and widespread swamp meadows.

4. Conclusions

In this paper, we took the cold and high-altitude SNP as a typical research area, the land-use/vegetation types of the three parks in SNP were extracted by multi-temporal Landsat-8 OLI multispectral remote sensing images, and cloud-free remote sensing images were obtained by using the de-clouding/snow/glacier algorithms. The verification sampling points were trained with a SVM classifier by combining field GPS sampling points and UAV sample data. Eventually, the alpine grasslands in cold and high altitudes were classified, and the data accuracy and validation were examined in this paper. We obtained the following conclusions:

(1)

The classification results showed that the SVM classifier could better distinguish the four vegetation types of alpine grassland, and the overall classification accuracy was relatively high. The average Kappa coefficient was 0.8366 among the vegetation types of the three parks, and the overall classification mapping accuracy could reach 84.52%. The mapping accuracy for the user was as high as 87.67%. The classification performances between the major land-use types were relatively high and the classification accuracies were very high. However, in the alpine grassland subcategories, the classification accuracies of the four typical grasslands were relatively low, especially between desert steppe and alpine meadow, and desert steppe and alpine steppe. This is highly related to the limited number of GPS sample points in the field, the low heterogeneity of the selected sample points, and high consistencies of the spectral characteristics of the four grassland types. The relatively low classification accuracy indicates the limitations of Landsat-8 multispectral remote sensing imageries in finer-resolution grassland classifications of high-altitude alpine mountains;

(2)

The accuracy comparison between the dataset produced in this paper and the other similar products demonstrated that the classification results of this paper were good at distinguishing the four typical types of alpine grasslands. In contrast, the other similar products tested were unable to distinguish grassland types either due to the low spatial resolution or inadequate classification system. This demonstrates that the datasets produced in this paper had the advantages of finer alpine grassland types and higher spatial/temporal resolutions;

(3)

The method in this paper can be applied to other similar cold and high altitudes with short vegetation growing seasons, but abundant clouds and snow/glaciers. The method in this paper can improve the efficiency of producing grassland type datasets and engineeringly generate year-by-year alpine grassland cover datasets, and provide high-quality data with high time efficiency and high spatiotemporal resolutions. The method can be utilized for other multispectral satellite imageries with the same band matching, such as Landsat 7, Landsat 9, Sentinel-2, etc. The results of this paper can facilitate further essential research on alpine grassland types, distribution, above-ground biomass, carrying capacity, and grassland degradation on the QTP with finer spatial and temporal resolutions.