Testing Accuracy of Land Cover Classiﬁcation Algorithms in the Qilian Mountains Based on GEE Cloud Platform

: The Qilian Mountains (QLM) are an important ecological barrier in western China. High-precision land cover data products are the basic data for accurately detecting and evaluating the ecological service functions of the QLM. In order to study the land cover in the QLM and performance of different remote sensing classiﬁcation algorithms for land cover mapping based on the Google Earth Engine (GEE) cloud platform, the higher spatial resolution remote sensing images of Sentinel-1 and Sentinel-2; digital elevation data; and three remote sensing classiﬁcation algorithms, including the support vector machine (SVM), the classiﬁcation regression tree (CART), and the random forest (RF) algorithms, were used to perform supervised classiﬁcation of Sentinel-2 images of the QLM. Furthermore, the results obtained from the classiﬁcation process were compared and analyzed by using different remote sensing classiﬁcation algorithms and feature-variable combinations. The results indicated that: (1) the accuracy of the classiﬁcation results acquired by using different remote sensing classiﬁcation algorithms were different, and the RF had the highest classiﬁcation accuracy, followed by the CART and the SVM; (2) the different feature variable combinations had different effects on the overall accuracy (OA) of the classiﬁcation results and the performance of the identiﬁcation and classiﬁcation of the different land cover types; and (3) compared with the existing land cover products for the QLM, the land cover maps obtained in this study had a higher spatial resolution and overall accuracy.


Introduction
Due to a combined influence of natural processes and human activities, the global land cover is changing rapidly [1]. With the advancement in multisource remote sensing data, land cover and its changes have been closely studied in recent years [2]. Accurate and reliable land cover data is indispensable, as it provides basic information related to scientific research in areas such as agricultural, environmental protection, and global change [3][4][5][6]. However, accurate land cover mapping and its changes are still facing many challenges due to land surface heterogeneity and spectral confusion, especially in higher-resolution mapping.
Since the late 1990s, many land cover data products with different resolutions have been developed, including the University of Maryland (UMD) classification [7], the International Global Biosphere Programme (IGBP) DISCover [8], Global Land Cover 2000 (GLC2000) [9], MCD12Q1 [10], GlobCover [11], and the Climate Change Initiative-Land Cover (CCI-LC) [12]; however, the resolution of these land cover products are low, ranging the different remote sensing classification algorithms in the QLM. Therefore, it is necessary to use different remote sensing classification algorithms to study the QLM, produce more accurate high-resolution land cover products, and assess the performances of the different classification algorithms.
According to the research background presented above, and based on the GEE cloud platform, Sentinel-2 remote sensing images were used to produce the land cover classification products for the QLM, with a spatial resolution of 10 m. Specifically, the research content of this study is as follows: (1) to evaluate the performance of the GEE cloud platform in land cover research in the QLM; (2) to assess the performances and accuracies of the different remote sensing classification algorithms in the land cover classification of the QLM; (3) to analyze the impacts of different characteristic variable combinations on the classification results while participating in classification; and (4) to compare and analyze the classification results obtained in this study using existing land cover products.

Study Area
The Qilian Mountains (35 • 48 -40 • 05 N, 93 • 18 -103 • 54 E) are located in northwestern China, to the northeast of the Qinghai-Tibetan Plateau (QTP), in Gansu and Qinghai Provinces ( Figure 1). The total area is 1.84 × 10 5 km 2 , with an elevation of 1623-5766 m [30]. The QLM is an important ecological security barrier in western China. Many inland rivers that originate from the QLM are important freshwater sources for the Hexi Corridor, maintaining the freshwater balance and oasis stability in the area. The QLM are composed of a number of NW-trending high mountains and valleys. The overall terrain characteristic is that the west is higher than the east, with high altitudes and complex and special geomorphological features. The QLM have a typical continental plateau climate, which is dry and cold, and it gradually becomes colder and drier from east to west. The annual average temperature is low, about 0.6 • C; and the ranges of annual and daily temperature are relatively large. The annual precipitation is about 400-700 mm, and it decreases from east to west and increases with elevation. The region has a high solar radiation intensity, with most areas receiving more than 2800 h of sunshine. The land cover types in the QLM are diverse, mainly including croplands, forests, grasslands, shrublands, wetlands, water bodies, construction lands, and bare lands, among which the grasslands and meadows account for a large percentage of the total area. The main crops are wheat, oilseed rape, highland barley, oats, maize, broad beans, and potatoes. classification results of the QLM, lacking a performance comparison of the different remote sensing classification algorithms in the QLM. Therefore, it is necessary to use different remote sensing classification algorithms to study the QLM, produce more accurate high-resolution land cover products, and assess the performances of the different classification algorithms. According to the research background presented above, and based on the GEE cloud platform, Sentinel-2 remote sensing images were used to produce the land cover classification products for the QLM, with a spatial resolution of 10 m. Specifically, the research content of this study is as follows: (1) to evaluate the performance of the GEE cloud platform in land cover research in the QLM; (2) to assess the performances and accuracies of the different remote sensing classification algorithms in the land cover classification of the QLM; (3) to analyze the impacts of different characteristic variable combinations on the classification results while participating in classification; and (4) to compare and analyze the classification results obtained in this study using existing land cover products.

Study Area
The Qilian Mountains (35°48′-40°05′N, 93°18′-103°54′E) are located in northwestern China, to the northeast of the Qinghai-Tibetan Plateau (QTP), in Gansu and Qinghai Provinces ( Figure 1). The total area is 1.84 × 10 5 km 2 , with an elevation of 1623-5766 m [30]. The QLM is an important ecological security barrier in western China. Many inland rivers that originate from the QLM are important freshwater sources for the Hexi Corridor, maintaining the freshwater balance and oasis stability in the area. The QLM are composed of a number of NW-trending high mountains and valleys. The overall terrain characteristic is that the west is higher than the east, with high altitudes and complex and special geomorphological features. The QLM have a typical continental plateau climate, which is dry and cold, and it gradually becomes colder and drier from east to west. The annual average temperature is low, about 0.6 °C; and the ranges of annual and daily temperature are relatively large. The annual precipitation is about 400-700 mm, and it decreases from east to west and increases with elevation. The region has a high solar radiation intensity, with most areas receiving more than 2800 h of sunshine. The land cover types in the QLM are diverse, mainly including croplands, forests, grasslands, shrublands, wetlands, water bodies, construction lands, and bare lands, among which the grasslands and meadows account for a large percentage of the total area. The main crops are wheat, oilseed rape, highland barley, oats, maize, broad beans, and potatoes.   The MultiSpectral Instrument (MSI) equipped on the Sentinel-2 satellite is a highresolution multispectral imaging sensor. The Sentinel-2 mission comprises a constellation of two polar-orbiting satellites (2A and 2B). Sentinel-2 has an inter-resolution of 5d after the two satellites are networked, and carries a multispectral imager. The spatial resolution of the visible light bands and near-infrared band is 10 m [31]. The remote sensing data used in this study were the Sentinel-2 Level-2A products for the 2020 and 2021 plant growth season (from the beginning of June to the end of August). A total of 899 Sentinel-2 images were selected in QLM, and each image included 12 spectral bands and 3 quality assessment (QA) bands. The data acquisition and preprocessing were conducted through online code writing on the GEE cloud platform. The QA60 band of the images was used and set as a cloud mask to remove the clouds from the images in order to obtain a cloud-free image.

Sentinel-1 Image Data
The Sentinel-1 data are C-band synthetic aperture radar (SAR) data [32]. Sentinel-1 images provided by the GEE cloud platform have undergone preprocessing procedures such as thermal noise removal, radiation calibration, and terrain correction. A total of 398 Sentinel-1 images that were synchronized with Sentinel-2 data were used in this study, and the data acquisition and processing were conducted using the GEE platform.

SRTM Data
The Shuttle Radar Topography Mission (SRTM) data were jointly measured by the National Aeronautics and Space Administration (NASA) and the National Imagery and Mapping Agency (NIMA) [33]. The product name of the digital elevation data provided in the GEE cloud platform was SRTMGL1_003, and the data acquisition and processing were conducted using the GEE platform.

Land Cover Datasets
The GlobeLand30 is a global land cover product produced by the National Geomatics Center of China [16]. The FROM-GLC30 is the first global 30 m spatial resolution land cover product, produced by Tsinghua University in 2013 based on Landsat images [17]. The FROM-GLC10 data were the latest results for the 2017 land cover produced by Tsinghua University, which used Sentinel-2 images and previous training samples based on Landsat data [18]. The Qilian Mountains 30 m land cover classification product data set (1985-2019) V2.0 (LCD-QLM (V2.0)) was also used in this study as auxiliary reference data [28]. The data were downloaded from the National Qinghai-Tibetan Plateau Science Data Center ( http://data.tpdc.ac.cn, accessed on 24 September 2021). The basic information about the land cover products is presented in Table 1.

Methods
Based on the GEE cloud platform and the remote sensing classification algorithms and remote sensing data it provided, this research aimed to classify the land cover in the Qilian Mountains. The research methods included the collection and optimization of the sample data, the construction of the feature space, and the machine learning classification algorithms and evaluation of classification results. The data processing and analysis flow diagram of this research is shown in Figure 2.

Methods
Based on the GEE cloud platform and the remote sensing classification algorithms and remote sensing data it provided, this research aimed to classify the land cover in the Qilian Mountains. The research methods included the collection and optimization of the sample data, the construction of the feature space, and the machine learning classification algorithms and evaluation of classification results. The data processing and analysis flow diagram of this research is shown in Figure 2.

Sampling Strategies
The collection of accurate training and validation samples is a necessary condition for accurate land cover classification [34]. Unrepresentative and/or inadequate samples will result in significant uncertainty of the land cover classification results [35]. Taking into account the actual situation in the QLM and referring to the existing land cover classification system and products, the classification system of this study is shown in Table 2.
The classification system used in this study was the land cover classification system based on ecosystem types, and was divided into 9 categories: croplands (CO), forests (FO), grasslands (GL), shrublands (SL), wetlands (WL), water bodies (WB), construction land (CL), bare land (BL), and permanent snow and ice (PSI). The high-resolution satellite imagery

Sampling Strategies
The collection of accurate training and validation samples is a necessary condition for accurate land cover classification [34]. Unrepresentative and/or inadequate samples will result in significant uncertainty of the land cover classification results [35]. Taking into account the actual situation in the QLM and referring to the existing land cover classification system and products, the classification system of this study is shown in Table 2. The classification system used in this study was the land cover classification system based on ecosystem types, and was divided into 9 categories: croplands (CO), forests (FO), grasslands (GL), shrublands (SL), wetlands (WL), water bodies (WB), construction land (CL), bare land (BL), and permanent snow and ice (PSI). The high-resolution satellite Remote Sens. 2021, 13, 5064 6 of 25 imagery of the Google Earth Engine was used as a base map; the sample points were selected from the areas where land cover changes have not occurred for many years, and the areas that were relatively uniform and had little interference from human activities to ensure the accuracy and authenticity of the sample data. Furthermore, the collected sample data were compared and verified with the 10 m ecosystem-type data of the Qilian Mountain Nature Reserve (part of the Qilian Mountains), which was obtained by artificial interpretation and field verification. The unit of each sample point was a pixel. The quality and quantity control of sample point data was realized by the method of area estimation adjustment error [36][37][38][39]. In this way, the deviation caused by stratified sampling could be adjusted. Using sample data, the confusion matrix n ij of the classification result could be calculated. A more informative presentation of the error matrix is the unbiased estimation p ij of the area ratio in the unit i and j of the error matrix: where W i is the area ratio of land cover types to the total area of the study area in the classification results, and the proportion of the area classified as type i is: where A tot is the total area of the map, and A i is the mapped area of land cover type i. An unbiased estimator of the total area of type j is then: where p j can be calculated using the following formula: The estimated standard error of the estimated area proportion is: The standard error of the error-adjusted estimated area is: An approximate 95% confidence interval for A j is: where the range of error is defined as the z-score multiplied by the standard error, and the value of the z-score is related to the confidence level (for 95% confidence, z = 1.96).
When the area of each land cover types in the classification results was within the estimated area ranges, it could be considered that the sample data on which the classification results were based were reasonable. The method accurately quantified the classification error caused by inaccurate sampling, and the final number of sample points is shown in Table A1 (Appendix A). Using the 'randomColumn' function in the GEE cloud platform, 70% of the sample data could be randomly selected for the classifiers training and image classification, and 30% of the sample data were used for the verification and evaluation of classification results.

• Spectral indices
Existing studies have shown that the application of remote sensing spectral indices can effectively improve the accuracy of the identification of different land cover types. In this study, the normalized difference vegetation index (NDVI) [40], the enhanced vegetation index (EVI) [41], the land surface water index (LSWI) [42], the normalized difference water index (NDWI) [43], the soil adjusted vegetation index (SAVI) [44], and the bare soil index (BSI) [45] were calculated using the GEE cloud platform based on the formulas for each spectral index, and each index was added to the original remote sensing image in turn (as shown in Figure 2). The calculation formulas for each index are as follows: Remote Sens. 2021, 13, 5064 where ρ SWIR , ρ NIR , ρ RED , ρ GREEN , and ρ BLUE are the surface reflectance values of the shortwave infrared, near-infrared, red, green, and blue bands, respectively. L is the soil regulation factor, which ranges from 0 to 1 and is usually assigned a value of 0.5 to better reduce the background difference of the soil and eliminate the noise impact of the soil [44].

Texture features
The texture features are an important attribute of remote sensing images, and different land cover types have different texture features. Based on the texture features, the accuracy of the recognition and classification can be improved [46]. The Gray-Level Co-occurrence Matrix (GLCM) is a classic method to extract the texture features [47]. It extracts textures through the conditional probability density between gray levels of remote sensing images, and is diffusely used in land cover classification research. The calculation of the GLCM can be obtained through the "glcmTexture" function in the GEE cloud platform. In this study, several commonly used texture features were selected for use in the classification, including the feature parameters of angular second moment (asm), contrast (con), correlation (corr), variance (var), inverse different moment (idm), sum average (savg), and entropy (ent) based on NDVI. The calculation of each feature variable was based on GLCM, they were dimensionless, and their range of values was not completely uniform.

Radar features
Studies have shown that SAR data are sensitive to land cover types such as water bodies, construction lands, and croplands [23]. The radar variables involved in the construction of the feature variables in this study included the backscatter coefficients of the Sentinel-1 data in VV polarization and VH polarization modes.

•
Terrain features In the recognition and classification of remote sensing images, the participation of terrain features can improve the accuracy of the classification. Therefore, based on the SRTMGL1_003 digital elevation data product, the "ee.Terrain.products" function provided by GEE was used to calculate the aspect, hill shade, slope, and elevation, and they were added to the remote sensing images as independent features.
According to the feature variables described above, the classifications of image were conducted in five feature variable combinations. In the first input feature variable combination, only spectral bands of Sentinel-2 participated in the classification of the image, and then the feature variable combinations of spectral indices, terrain features, radar features, and texture features were added in turn.

Classification Algorithms
At present, remote sensing classification methods such as the support vector machine (SVM), the classification regression trees (CART), and the random forest (RF) algorithms have been widely used in land cover mapping and crop type identification [4]. In order to analyze the differences in the accuracies of the extraction and classification of the different land cover types when using different classification algorithms, and when using the same feature variable combination and the differences in the accuracies of the extraction and classification of different land cover types when using the same classification algorithm but different combinations of feature variables, three classification algorithms were used in this study. These three classification algorithms are introduced below.

• Support Vector Machine
The support vector machine (SVM) was proposed by Vapnik in 1995 [48]. The SVM has significant advantages in terms of nonlinear problems, small samples, and a high dimensionality. It is widely used due to its small training samples and support for highdimensional feature spaces. The parameters that need to be adjusted when using the SVM to classify remote sensing images are the kernel function type, the gamma value of the kernel function, and the cost parameter.

• Classification and Regression Tree
The classification and regression tree (CART) was proposed by Breiman et al. in 1984 [48]. It is widely used in land cover extraction and remote sensing image classification research due to its simple structure, fast calculation speed, and the advantage of being easy to understand. The parameters that need to be optimized when using CART to classify remote sensing images are the maximum and minimum numbers of leaf nodes in each tree.

•
Random Forest The random forest (RF) algorithm was proposed by Leo Beriman in 2001 [49]. Studies have shown that the RF has the advantages of stability, rapidity, and high precision in processing remote sensing data; and therefore, it has become the most widely used classifier among current remote sensing classification algorithms. It has important applications in crop extraction, image classification, and agricultural regression models [50,51]. The two main parameters that need to be adjusted and optimized for the RF when classifying remote sensing images are the number of decision trees that need to be created and the minimum number of leaf nodes. Studies have indicated that the values of the RF parameters have little effect on the accuracy [51], so the number of decision trees was set as 300. Additionally, when using RF for classification in GEE, the importance score of the parameters participating in the classifier can be calculated by using the "explain" function. The value of the importance score is not absolute, and there is no uniform and fixed value range, but changes with the number of sampling data and feature parameters participating in the classification [52,53]. The importance and contribution of the feature parameters participating in the classification can be determined by the relative value of the importance score in a specific situation.

Accuracy Assessment
The confusion matrix is a method that is commonly used to assess the accuracy of image classifications [54]. In this study, the calculation of the confusion matrix was conducted through online programming of the GEE cloud platform, and then the overall accuracy (OA), Kappa coefficient, producer accuracy (PA), and user accuracy (UA) were calculated. Among them, the OA and the Kappa coefficient could fully reflect the comprehensive accuracy of the results, and the PA and UA could be used to assess the classification accuracy of a specific land cover type. The calculation formulas were as follows: where N is the total number of samples used for the accuracy assessment; n is the total number of columns in the confusion matrix; X ii is the number of samples in the i-th row and i-th column in the confusion matrix; and X i+ and X +i are the total number of samples in the i-th row and i-th column, respectively. The uncertainty analysis of the accuracy of classification results is realized by calculating the error range, which is defined as the z-score multiplied by the standard error [36].
In this study, all of the data processing and calculations, including data acquisition, data processing, composition of images, construction of feature space, calculation of parameters, implementation of classifiers, and calculations of confusion matrices, were all implemented by the GEE JavaScript API.

Classification Results and Accuracy of Classification Results
Remote sensing classification algorithms, including SVM, CART, and RF, were used to classify the Sentinel-2 composite images in the QLM and to analyze the ability of the different classifiers to identify land cover types in high spatial resolution images. The land cover results based on the three classification algorithms are shown in Figure 3.
In this study, all of the data processing and calculations, including data acquisition, data processing, composition of images, construction of feature space, calculation of parameters, implementation of classifiers, and calculations of confusion matrices, were all implemented by the GEE JavaScript API.

Classification Results and Accuracy of Classification Results
Remote sensing classification algorithms, including SVM, CART, and RF, were used to classify the Sentinel-2 composite images in the QLM and to analyze the ability of the different classifiers to identify land cover types in high spatial resolution images. The land cover results based on the three classification algorithms are shown in Figure 3.  The results showed that the classification results obtained using the different remote sensing classification algorithms had high consistency in terms of the proportions and the distribution characteristics of the land cover types. Among them, bare land and grassland accounted for the largest proportion of the study area, accounting for more than 90%.
The grassland was mainly distributed in the eastern and central regions of the QLM, and the relatively warm and humid valleys in the west were also distributed within a small area. The bare land was widely distributed in the western part of the QLM, where the temperature is low and the precipitation is scarce. In addition, the cropland was concentrated in the eastern basin, the middle reaches of the Datong River, and the plains on the northern slope of the northeastern part of the mountains, where the land is flatter, with abundant water resources and suitable temperatures. The forests were distributed in the valleys in the middle and eastern parts of the study area and on the northern slopes of the QLM in the northeast. The permanent snow and ice were distributed in the mountainous areas with higher elevations (4891 m on average) in the western and central parts of the study area.
The confusion matrix of the classification results obtained using the different remote sensing classification algorithms is shown in Table 3. In general, in multiple experiments, the OAs of the three classifiers reached 95%, which indicated that the three classification algorithms were suitable for high-spatial resolution image classification in the study area. Among them, the average value of OAs of the RF was 96.51% and the Kappa coefficient was 0.95, and the OAs of the SVM and CART were 94.67% and 94.50% and the Kappa coefficients were 0.92 and 0.91, respectively. This indicated that the OA of the classification results generated by the RF was higher than those of the classification results of the CART and SVM, which showed that the RF classifiers could more accurately identify and classify the land cover types in the remote sensing images. Specifically, all three classification algorithms had high classification accuracies for grassland, water bodies, bare land, and permanent snow and ice, but their classification accuracies for cropland, forests, and construction land were significantly different. The PA and OA of the RF for cropland and construction land were higher than those of the other two classifiers. Compared with the CART and RF, the SVM had a lower recognition and classification accuracy for forests. Using the "explain" function provided by the GEE cloud platform, the importance of each of the feature variables of the random forest (RF) classifier when participating in the classification was determined. The contribution of the variables to the classification results was greater if they had a higher importance score [55]. Figure 4 shows the importance score distribution of the 31 input feature variables involved in the classification. As can be seen from the figure, terrain features and radar features had higher importance scores, while the lower importance scores were for spectral bands, spectral indices, and texture features. Specifically, among the terrain features, the importance of the elevation feature was the highest (801.24), followed by slope feature (789.51), and the importance of the aspect and hill shade were relatively low. Radar features including VV and VH had high importance scores, with average values of 769.23 and 681.5, respectively. Spectral bands B1, B11, B5, B12, B3, and B2 of the Sentinel-2 images had generally higher importance scores during the classification, and their average values were all >560. Among the six spectral indices, the LSWI had the highest importance scores, with an average value of 603.5, which indicated that they were very important for water recognition. The importance scores of the texture features were generally low, among which the NDVI_idm, NDVI_corr, NDVI_asm, and the NDVI_ent were the lowest, with their values all <560. = Kappa coefficient.

Influence of the Feature Variables on the Classification Accuracy
3.2.1. Importance Scores of the Variables Used in the RF Classification Algorithm Using the "explain" function provided by the GEE cloud platform, the importance each of the feature variables of the random forest (RF) classifier when participating in t classification was determined. The contribution of the variables to the classification resu was greater if they had a higher importance score [55]. Figure 4 shows the importan score distribution of the 31 input feature variables involved in the classification. As can seen from the figure, terrain features and radar features had higher importance scor while the lower importance scores were for spectral bands, spectral indices, and textu features. Specifically, among the terrain features, the importance of the elevation featu was the highest (801.24), followed by slope feature (789.51), and the importance of t aspect and hill shade were relatively low. Radar features including VV and VH had hi importance scores, with average values of 769.23 and 681.5, respectively. Spectral ban B1, B11, B5, B12, B3, and B2 of the Sentinel-2 images had generally higher importan scores during the classification, and their average values were all >560. Among the spectral indices, the LSWI had the highest importance scores, with an average value 603.5, which indicated that they were very important for water recognition. The i portance scores of the texture features were generally low, among which the NDVI_id NDVI_corr, NDVI_asm, and the NDVI_ent were the lowest, with their values all <560.

Influence of the Feature Variables on the OA
A comparison of the overall classification accuracies of the different classifiers when different feature variables were inputted is shown in Figure 5. It can be seen that the OA of the same classifier was different when different feature variables were involved in the classification. For the SVM, in addition to the texture features, the participation of other feature variables improved the value of OA. In particular, the participations of spectral indices and radar features significantly improved the OA of SVM classifier, with a value of 0.76% and 1.33%, respectively. The SVM classifier had the highest OA when the spectral bands, spectral indices, terrain features, and radar features participated in the classification at the same time, with an OA of 95.65% and a Kappa coefficient of 0.93. The participations of each type of feature variables increased the OA of the classification results when CART and RF were used for classifications, respectively. They both reached the highest value of OAs when the feature variables were combined with the spectral bands, spectral indices, terrain features, radar features, and texture features. The highest OA of the CART was tion at the same time, with an OA of 95.65% and a Kappa coefficient of 0.93. The parti pations of each type of feature variables increased the OA of the classification results wh CART and RF were used for classifications, respectively. They both reached the highe value of OAs when the feature variables were combined with the spectral bands, spect indices, terrain features, radar features, and texture features. The highest OA of the CAR was 95.44%, and the Kappa coefficient was 0.93. For the RF, the highest OA was 97.18 and the Kappa coefficient was 0.95. Figure 5. The influences of the different feature variable combinations on the overall accuracy (OA) of the classification results. Feature variable combination 01: spectral bands; feature variable combination 02: spectral bands + spectral indices; feature variable combination 03: spectral bands + spectral indices + terrain features; feature variable combination 04: spectral bands + spectral indices + terrain features + radar features); feature variable combination 05: spectral bands + spectral indices + terrain features + radar features + texture features.

Influence of the Feature Variables on the PA of the Different Land Cover Types
Table A2 (Appendix A) shows the producer accuracies (PAs) of the different lan cover types obtained using the different remote sensing classification algorithms wh different feature variables participated in the classification. The PAs of the forests, gra lands, water bodies, bare lands, and permanent snow and ice were generally higher a were less affected by the inputted feature variable combinations, while the land cov types of croplands, shrublands, wetlands, and construction land were greatly affected the feature variable combinations. Specifically, when the RF classifier was used and t feature variable combination 03: spectral bands + spectral indices + terrain features; feature variable combination 04: spectral bands + spectral indices + terrain features + radar features); feature variable combination 05: spectral bands + spectral indices + terrain features + radar features + texture features.

Influence of the Feature Variables on the PA of the Different Land Cover Types
Table A2 (Appendix A) shows the producer accuracies (PAs) of the different land cover types obtained using the different remote sensing classification algorithms when different feature variables participated in the classification. The PAs of the forests, grasslands, water bodies, bare lands, and permanent snow and ice were generally higher and were less affected by the inputted feature variable combinations, while the land cover types of croplands, shrublands, wetlands, and construction land were greatly affected by the feature variable combinations. Specifically, when the RF classifier was used and the feature variable combination was spectral bands + spectral indices + terrain features + radar features, the PA value of the cropland reached the highest value (71.25%). When all of the feature variables were inputted into the RF classifier, the PA of the forest was the highest (94.13%). The PA of the grassland reached the maximum value of 98.61% when only the spectral bands were inputted into the SVM classifier. When the CART classifier was used and the feature variable combinations was spectral bands + spectral indices + terrain features, the PAs of the shrublands and water bodies reached the highest values (66.67% and 97.03). When the RF classifier was used and the spectral bands, spectral indices, terrain features, and radar features were involved in the classification, the PA of the wetlands was the highest (49%). When the SVM classifier was used and all of the feature variables were involved in the classification, the PA of the construction land reached the highest value of 90%. When the spectral bands and spectral indices were inputted into the SVM classifier at the same time, the PA of the bare lands reached the highest value of 99.4%. The PA of the permanent snow and ice reached the maximum value of 100% when only the spectral bands were inputted into the SVM classifier.

Influence of the Feature Variables on the UA of the Different Land Cover Types
Table A3 (Appendix A) shows the user accuracies (UAs) of the different land cover types for different remote sensing classification algorithms when different feature variables participated in the classification. The UA of the croplands reached the maximum (100%) when only the spectral bands were inputted into the SVM classifier. When the RF classifier was used and all of the feature variables were involved in the classification, the UA of the forest reached the highest value (94.15%). The UA of the grassland reached the maximum value (96.97%) when the spectral bands, spectral indices, terrain features, and radar features were inputted into the RF classifier. There were eight cases in which the UA of shrublands reached the highest value of 100%. When the SVM classifier was used and only the spectral bands were involved in the classification, the UAs of the wetlands, water bodies, and construction lands reached their highest values (100%, 99.21%, and 100%, respectively). The UA of the bare land reached the maximum value (98.63%) when all of the feature variables were inputted into the RF classifier. When the SVM classifier was used and the spectral bands, spectral indices, and terrain features were involved in the classification, the UA of the permanent snow and ice cover reached the highest value (99.79%).

Comparison of Classification Results with Other Land Cover Products
In order to assess the land cover results obtained in this study, the classification results obtained in this research were compared with existing high-spatial-resolution land cover products, including the FROM-GLC30, FROM-GLC10, Land Cover Dataset for the Qilian Mountains Area from 1985 to 2019 (V2.0) (LCD-QLM (V2.0)), and GlobeLand30. Figure 6 presents a visual comparison of the classification results obtained in this research and the existing land cover products in three different magnified areas. The biggest difference between the different land cover products lay in areas with complex land cover types. Therefore, three areas in the QLM with complex land cover types that included construction land, grassland, forests, farmland, and water bodies were selected as the magnified areas. We also took the standard false-color composited images of the three regions as the true value as a reference (Figure 6a-c). It can be seen that the classification effects of the land cover results obtained in this research were significantly improved. For example, in term of the identification and classification of construction lands, the zoomed-in display of the classification results obtained in this study (Figure 6d-f) were more sufficient than FROM-GLC10 (Figure 6j-l), although FROM-GLC10 had a higher spatial resolution and accuracy among many land cover products. In addition, compared with the FROM-GLC30, GlobeLand30, and LCD-QLM (V2.0) products, the classification results of this research had a higher spatial resolution and more accurate identification of the boundaries between the land cover types.
For the proportions of the total QLM area occupied by the different land cover types, several land cover products were used to analyze the consistency of and differences between the different land cover products (Table 4). In general, the distribution of the land cover types for several of the land cover products was basically consistent, mainly for bare land and grassland, and the sum of their area accounted for about 88% of the total area of the QLM. However, specifically, the proportions of some land cover types in the study area were different. Among them, the differences in the proportions of the grassland and unused land were more obvious, and the differences in the proportions of the water bodies and construction land were relatively small. Compared with other products, the proportion of the cropland in LCD-QLM (V2.0) was significantly smaller. The proportions of the grassland and permanent snow and ice in the study area decreased with increasing spatial resolution, while the proportions of bare land increased with increasing spatial resolution. this research had a higher spatial resolution and more accurate identification of the boundaries between the land cover types.  The same validation sample dataset was used to evaluate the accuracies of all of the land cover products to analyze whether the results obtained in this research were better than the existing high-spatial-resolution global land cover products. Table 5 compares the accuracies of the classification results obtained in this study with those of the existing land cover products. It can be seen that the overall accuracies and Kappa coefficients of the products obtained in this study were higher than those of the data products of the FROM-GLC30, FROM-GLC10, LCD-QLM (V2.0), and GlobeLand30. Therefore, the classification effect of the 10 m spatial resolution land cover maps of the QLM obtained in this study was better than those of the existing high-resolution land cover products. Table 5. Accuracy comparison (%) of the classification results obtained in this study with the existing land cover products.

Comparison of the Performances of the Different Classification Algorithms
Comparing the performances of the different classification algorithms used for the land cover classification of QLM was one of the important research goals of this study. Three classification algorithms, including the support vector machine (SVM), classification regression tree (CART), and random forest (RF) algorithms were selected for analysis in this research, and the classifications were performed using the feature variable combinations with the best classification effect. The classification accuracies were compared, and it was found that the three remote sensing classification algorithms produced accurate classification results. Among them, the RF had the best OA, followed by the SVM and the CART. This result was inconsistent with the results of several existing studies. For example, Abdi used Sentinel-2 remote sensing data and remote sensing classifiers, including the SVM, RF, extreme gradient boosting (Xgboost), and deep neural network (DNN) algorithms, for land cover classification mapping of a 10 km × 12 km area in Uppsala County, south-central Sweden, and determined that the SVM had the highest OA, followed by the Xgboost and RF, and the DNN had the lowest OA [56]. Rana et al. used Sentinel-2 remote sensing data to compare the effectiveness of the traditional and principal component analysis (PCA)based methods of different classification algorithms, including the maximum likelihood estimation (MLE), RF, and SVM. Their results showed that whether the PCA was used, the SVM had a higher accuracy than the RF [57]. These results were also similar to the previous reported studies. For example, Talukdar et al. used Landsat 8 remote sensing data and classification methods such as the artificial neural network (ANN), SVM, fuzzy ARTMAP (FA), RF, Mahalanobis distance (MD), and spectral angle mapper (SAM) for classification, and their results showed that the RF had the highest classification accuracy, followed by the ANN, SVM, fuzzy ARTMAP, and SAM; and the MD was the lowest [4]. Ge et al. used Landsat 8 image data and remote sensing classification algorithms, including the k-nearest neighbor (KNN), RF, SVM, and ANN, to carry out land cover classification research in the Dengkou Oasis in China. They found that the ANN had the highest classification accuracy, followed by the RF, SVM, and KNN [58]. The basic reasons for the difference in results were the different study area, topographic characteristics, and inconsistencies in the land cover types. Secondly, different satellite data sources were used in the different studies, and even the same data source with different time phases can also cause changes in the classification results. Therefore, the differences in the study areas and remote sensing data sources participating in the classification had a greater influence on the effects of the different classifiers, and the classification effect of the same remote sensing classification algorithm varied from one region to other.

Influence of Feature Variables on Remote Sensing Classification
The GEE cloud platform can be used to evaluate the contribution of the feature variables participating in the RF classifier. In this study, the average importance scores from high to low were as follows: terrain features, radar features, spectral bands, spectral indices, and texture features. This result was inconsistent with the results of several existing studies. For example, Li et al. used the RF classifier and fused Landsat data and Sentinel data to classify and map Africa. Their results indicated that the highest importance scores for the feature variables participating in the classification were those of the NDVI_S2 and NDVI_L8, and the lowest importance scores were for Band2 and Band8A of the Sentinel-2 and Band2 and Band10 of Landsat 8 [52]. Zhang et al. used Landsat 8 OLI data and a GEE-based RF classifier to map the urban areas of three cities in China. In the importance analysis of the variables, NDVI and NDWI had the highest importance scores [59]. These results were also consistent with the findings of previous reported studies. Liu et al. used the RF classifier to classify the Landsat TM and OLI data for the Gannan Prefecture from 2000 to 2018. Their results showed that the topographic features, including the altitude and slope, contributed the most, followed by the spectral indices and spectral bands [53]. Phan et al. used Landsat 8 data from different seasons and the RF classifier based on GEE to classify the land cover in Mongolia. The results showed that among all the input variables, elevation was the most important variable, followed by other feature variables [60]. Therefore, if the data sources used for the classification and the study areas were different, the importance scores of the feature variables input into the RF in the different studies were also different.
When different classification algorithms were used for classification, the participations of different feature variables had different impacts on the PA values of the different land cover types (Table A2 (Appendix A)). When three classification algorithms were used in an image classification, compared with the land cover types of the croplands, shrublands, wetlands, and construction lands, other land cover types had high PAs and were less affected by the feature variable combinations. In addition, compared with CART and RF, the differences in the feature variables had a greater impact on the PA of the various land cover types when the SVM classification algorithm was used. Specifically, the participation of the most feature variables in the classification increased the PA values of the land cover types, including croplands, shrublands, wetlands, and construction lands. It is worth noting that the participation of spectral indices greatly increased the PA of croplands, by 24.91%; and the participations of terrain features and radar features greatly increased the PA of construction lands, by 45.02% and 41.09%, respectively. In contrast, there were also cases in which the feature variables had negative impacts on the PAs of some of the land cover types. Among them, the participation of texture features caused decreases in the PAs of the croplands and construction lands, by 3.27% and 1.66%, respectively. For the CART classification algorithm, in general, the participation of most of the feature variables had a positive effect on the increase in the PAs of most of the land cover types. Among them, the participation of the terrain features significantly increased the PAs of wetlands and construction lands, by 17.91% and 13.64%, respectively; and the participation of the radar features significantly increased the PA of construction lands, by 14.24%. However, when the radar features were inputted, the PAs of the croplands, shrublands, and wetlands decreased slightly (by 5% on an average). When the RF algorithm was used, the participations of the terrain features and radar features had a greatly positive effect on the PAs of the wetlands and construction lands. However, the participation of the spectral indices had a small negative effect on the PAs of some land cover types, including the shrublands, wetlands, and construction lands, with an average decrease of 2-7.83%.
When the different classification algorithms were used for the classification, the participation of the different feature variables also had different impacts on the UA (Table A3 (Appendix A)). When the SVM classification algorithm was used for the classification, the participation of most of the feature variables had a positive effect on the improvement of the UA values of the land cover types. Moreover, the most significant effect was that the UAs of the wetlands and construction lands increased by 8.88% and 19.86%, respectively, when the radar features were inputted into the classification process. In contrast, there were also cases in which the feature variables had a negative effect on the UAs of some of the land cover types. Specifically, the participations of the spectral indices and texture features resulted in a significant decrease in the UA of the wetlands (16.67% and 28.01%). When the terrain features participated in the classification process, the UA of the construction lands decreased by 27.8%. For the CART classification algorithm, in general, the participation of most of the feature variables had a positive effect on the improvement of the UAs of the land cover types. The most substantial effect was that the UAs of the croplands and shrublands, which increased by an average value of 12%, and the UA of wetlands increased by 23.02%. When the spectral indices, terrain features, and radar features were added to the classifier in turn, the UA of the construction lands gradually increased by 11.78%. However, when the spectral indices were added to the classifier, the UAs of the croplands and grasslands decreased by 8.56%. In addition, the participation of the radar decreased the UAs of the croplands, shrublands, and wetlands by an average value of 5.46%. For the RF classifier, the different land cover had higher UA values, and the feature variables had positive effects on the increase of the UA values, except for the wetlands. Specifically, when the terrain features, radar features, and texture features were added to the classifier in turn, the UA of the wetlands gradually decreased by 2.14%.

Comparison of the Land Cover Results Obtained in This Study with Existing Land Cover Products
In this research, specific land cover classifications were performed, and then compared with the existing global high-spatial-resolution products in the QLM. The land cover results obtained in this study had higher spatial resolutions and classification accuracies. Moreover, the identification and classification of the land cover types in a relatively small area and the boundaries between the land cover types were more accurate. Nevertheless, the classification results obtained in this study were different from the existing global high-spatial-resolution land cover products in terms of consistency. The data sources, classification algorithms, and feature variables involved in the classification were different, which were the main reasons for the differences in the classification results. Optical remote sensing imagery has the advantages of long-term coverage and wide coverage, and it is the main data source for land use classification research. Different data sources have different spatial and temporal resolutions, which leads to uncertainties in the classification results and inconsistencies between the different land cover products. In the existing high-spatial-resolution land cover products and land cover studies, most products used Landsat images with a resolution of 30 m as the data source, which was different from the Sentinel-2 images used in this study. For example, the GlobeLand30 developed by Chen et al. used Landsat TM/ETM+ data from 2000 and 2010, supplemented by HJ-1 satellite data [16]. For the FROM-GLC30 developed by Gong et al.,8929 Landsat TM/ETM+ scenes of the green season from 1984 to 2011 were collected from various sources [17]. Based on the Landsat 8 surface reflectance imagery from 2014 to 2016, Zhang et al. developed global land cover products, including GLC_FCS30-2015 [19]. Second, in addition to the supervised classification methods that are commonly used, some new methods have also been proposed and used in land cover classification. For example, Zhong et al. [61,62] proposed a new time series land cover mapping method based on machine learning and applied it to the development of land cover products in the Qilian Mountains. The Chinese land cover products developed by Zhang et al. were based on a new SPECLib-based operational approach proposed in their research [27]. The pixel-based, object-oriented, and knowledge fusion (POK)-based classification methods were used with the GlobeLand30 dataset [16]. Therefore, different classification algorithms will also result in differences in the land cover results.

Limitations and Prospects of Land Cover Classification in QLM
There are many factors that can lead to an error in land cover classification. The phenomenon of "the same subject with different spectra" and "different subjects with the same spectra" are the most vital reasons for the misclassification of land cover types. For example, misclassification of forest and grassland occurred in some areas in this study. Secondly, the limitations of the classifiers also lead to classification biases. Specifically, although the three methods achieved high classification accuracy, there were still some problems worth noting. For example, the CART classifier had the problem of confusion between bare land and water bodies. The SVM had a good recognition effect on urban areas, but it recognized a very small area of bare land as construction land. In general, the land cover types with the fewest misclassification phenomenon were water bodies and permanent snow and ice. Third, the problem of mixed pixels was also one of the sources of classification errors. In addition, the error source of classification was also caused by cloud removal and hill shadows. The cloud removal process was performed during the preprocessing of remote sensing images, but there were still very few cloud shadows and hill shadows in the composited image, which created the classification biases.
With the gradual enrichment of remote sensing data sources, the opportunities and challenges involved in obtaining land cover products using remote sensing image identification and classification are gradually increasing. The gradual improvement of the spatial resolution of satellite data will produce more land cover products with finer spatial resolutions. The application of multisource data fusion and the continuous update and improvement of the classification algorithms will greatly improve of the accuracies of the products. Therefore, remote sensing images with higher temporal and spatial resolutions can be used as data sources for land cover studies in the future, and the classification system can be further refined. In addition, shortening the time period of the products in the QLM, such as the production of land cover products on a quarterly or monthly basis, will be more conducive to monitoring the vegetation changes in the area, and will play a better auxiliary role in the studies of seasonal changes and mutation detection of vegetation in this area.

Conclusions
Based on the Google Earth Engine and the various remote sensing satellite data provided by it, a feature space composed of feature variables such as spectral bands, spectral indices, terrain features, texture features, and radar features was constructed. Three remote sensing image classification algorithms, including the SVM, CART, and RF algorithms, were used to automatically extract the land cover types in the Qilian Mountains. The accuracies of the classification results were evaluated, and the impacts of the different classification algorithms and the different feature variable combinations on the classification results were analyzed. In addition, the comparison between the classification results obtained in this research and existing high-spatial-resolution products were conducted. The results showed that the GEE platform, based on the cloud framework, had the advantages of convenient data acquisition and strong data-processing capabilities, and could be used to quickly and accurately extract and map land cover types in a large area. Moreover, the three remote sensing image classification algorithms assessed in this study could be used to obtain classification results with high classification accuracies. The classification results obtained using the three classification algorithms had high consistency in terms of the area ratio and the distribution characteristics of the land cover types, but there were still differences in the classification accuracies of the classification results. Among them, in terms of the OA of the classification results, the RF classifier was the highest, with an average OA of 96.51%, followed by the SVM classifier (94.67%) and CART classifier (94.50%). Therefore, the most appropriate classification algorithm was RF for the land cover classification in QLM. Furthermore, the participation of feature variables in the classification process improved the OA of the land cover classification results. Specifically, the participations of spectral indices and radar features significantly improved the OA of the SVM classifier, with values of 0.76% and 1.33%, respectively; the participations of spectral indices significantly improved the OA of the CART classifier, with a value of 1.65%; and the participations of each type of feature variables increased the OA of the RF classifier, with an average value of 0.33%. Nevertheless, there were still cases in which the different feature variables had different effects on the different classifiers, and the contributions of the different feature variables to the identification and extraction of the different land cover types were also different. These phenomena were reflected in the UA and PA values of different land cover types in the classification results when different feature variables combinations participated in the classification. Finally, compared with the existing global land cover products, the 10 m spatial resolution land cover map of the QLM obtained in this research had a higher spatial resolution and accuracy. When the same validation sample dataset was used to evaluate the accuracies of all land cover products, the OA of the classification results of this research was higher than the OA of FROM-GLC10, with a value of 7.51%. This gap was much higher compared to other land cover products. The results of this study provided a scientific basis for other related studies in the QLM, as well as a reference for high-spatial-resolution land cover classification and mapping of large areas.