Automatic Extraction of the Spatial Distribution of Picea schrenkiana in the Tianshan Mountains Based on Google Earth Engine and the Jeffries–Matusita Distance

Abstract

As a distinct species in the Tianshan Mountains (TS) of Central Asia (CA), Picea schrenkiana plays a significant role in water purification, soil and water conservation, and climate regulation. In the context of climate change, rapidly and accurately obtaining its spatial distribution has critical decision-making significance for maintaining ecological security in the arid area of CA and the sustainable development of the “Silk Road Economic Belt”. However, conventional methods are extremely challenging to accomplish the high-resolution mapping of Picea schrenkiana in the TS, which is characterized by a wide range (9.97 × 10⁵ km²) and complex terrain. The approach of geo-big data and cloud computing provides new opportunities to address this issue. Therefore, the purpose of this study is to propose an automatic extraction procedure for the spatial distribution of Picea schrenkiana based on Google Earth Engine and the Jeffries–Matusita (JM) distance, which considered three aspects: sample points, remote-sensing images, and classification features. The results showed that (1) after removing abnormal samples and selecting the summer image, the producer accuracy (PA) of Picea schrenkiana was improved by 2.95% and 0.24%–2.10%, respectively. (2) Both the separation obtained by the JM distance and the analysis results of eight schemes showed that spectral features and texture features played a key role in the mapping of Picea schrenkiana. (3) The JM distance can seize the classification features that are most conducive to the mapping of Picea schrenkiana, and effectively improve the classification accuracy. The PA and user accuracy of Picea schrenkiana were 96.74% and 96.96%, respectively. The overall accuracy was 91.93%, while the Kappa coefficient was 0.89. (4) The results show that Picea schrenkiana is concentrated in the middle TS and scattered in the remaining areas. In total, 85.7%, 66.4%, and 85.9% of Picea schrenkiana were distributed in the range of 1500–2700 m, 20–40°, and on shady slope and semi-shady slope, respectively. The automatic procedure adopted in this study provides a basis for the rapid and accurate mapping of the spatial distribution of coniferous forests in the complex terrain.

1. Introduction

As one of the most valuable terrestrial ecosystems, the forest ecosystem plays a vital role in maintaining global carbon balance and slowing down the rise in the greenhouse gas concentration [1]. The carbon stock of forest ecosystems accounts for about 33%–46% of terrestrial ecosystems [2]. Under the background of “carbon neutrality”, global changes caused by both natural and human factors directly or indirectly affect forest ecosystems, which has aroused widespread concern in the whole society [3]. The arid area of Central Asia (CA), which is connected to China by mountains and rivers, is a key hub zone in the construction of “the Silk Road Economic Belt” and a hot area concerned with the international community [4]. The Tianshan Mountains (TS) is one of the most important parts in the arid area of CA, in which the forest ecosystem has critical ecological values and functions such as carbon sink [5], water resource conservation, flood regulation, water purification, and soil and water conservation. The TS is highly essential for maintaining the ecological security of the downstream areas and the sustainable development of the whole arid zone [6]. The arid areas of CA, far from the sea and interspersed with mountains and basins, have light precipitation but intense evaporation, severe water shortages, and fragile ecosystems. Therefore, this region is a dramatic response region to global changes [7,8]. In the past half-century, although the trend of “warming and humidification” has occurred in the arid area of CA [9], it has not slowed down the ecological degradation in the region, because the positive ecological effect of increasing precipitation is not sufficient to offset the negative effect of rising temperatures [10]. In addition, the normalized difference vegetation index (NDVI) of the natural vegetation growing season in CA has declined at a rate of 0.0003 per year due to a sudden increase in the temperature since 1998 [10]. All of the above studies indicate that the natural vegetation in this area shows browning as a whole. Picea schrenkiana is a coniferous forest species with an altitude range of 1400–2800 m on the northern slopes of the TS, and is particularly sensitive to climate change [11]. Timely and accurate acquisition of its spatial distribution information is of great significance to ensure the sustainable development of the arid area of CA and “the Silk Road Economic Belt”.

Remote sensing is the most time-efficient and economical method for obtaining large spatial-scale land-cover information [12,13]. With the improvement in remote-sensing image resolution, remote-sensing technology plays a vital role in forest-cover extraction and forest-change monitoring [14]. In Vietnam, researchers utilized Landsat images from 1973 to 2020 to identify changes in the spatial distribution of mangroves in Thanh Hoa and Nghe An provinces, and analyzed the reasons for these changes [15]. The researchers first determined the range of cultivated land in the Guanzhong region based on the land-use map, and then distinguished between the grain/non-grain crops based on MODIS vegetation index data. Finally, the non-grain production of cultivated land was extracted [16]. The above studies on land-cover extraction and forest-change monitoring mainly used Landsat data with high spatial resolution and MODIS data with high temporal resolution. Landsat data are limited by the temporal resolution and rainy weather, and MODIS data are limited by the spatial resolution and mixed pixels, which when used separately do not meet the requirements for fine forest-cover extraction and forest-change monitoring. To solve this problem, some scholars [17,18] implement the method for integrating different remote-sensing data to make comprehensive use of the advantages of different sensors in time and space. However, the traditional data acquisition and data processing and analysis methods based on stand-alone off-line remote-sensing analysis software (ENVI and ERDAS) have great limitations [19], which cannot meet the needs of wide-area and high spatial resolution fine forest-cover extraction and forest-change monitoring.

The rapid development of remote-sensing cloud-computing platforms has become a new way to solve geographical and spatial-related problems. Researchers can directly invoke the huge amount of data on the platform, in conjunction with the algorithms provided by it, to implement online analysis and processing, as well as visualization operations, which considerably increase the efficiency of data analysis and processing. So far, Google Earth Engine (GEE) has been used in forestry [20], agriculture [21,22], wetlands [23], disaster monitoring [24,25], and land-use change [26,27]. In terms of vegetation information extraction, Liu et al. [28] selected Landsat images and used decision-tree classification to extract the spatial distribution information of mangroves in Yuxian County in Vietnam, and monitored the inter-annual dynamic changes in mangroves in the region from 1993 to 2017. Li et al. [29] selected Landsat OLI and MODIS NDVI data and fused spectral and phonological characteristics to extract distribution information on rubber forests in Thailand, Malaysia, and Indonesia. Xu et al. [30] used the GEE platform and multi-temporal Landsat images to extract information from citrus orchards in southern Jiangxi based on the random forest (RF) method, and the average OA reached 93.15%. The above results show that with the help of GEE, there are three main features for the extraction of vegetation information by domestic and foreign scholars. First, this vegetation is the most important cash crop in the region, such as rubber and citrus, which plays an essential role in the prosperity of farmers and the economic development of the region. Second, it has a major ecological value and ecological function. Third, it grows mainly in tropical and subtropical regions, such as Southeast Asia and Yunnan, Guangxi, Guangdong, and Fujian provinces.

At present, the research on Picea schrenkiana mainly includes nitrogen addition [31], individual genomics [32], stoichiometric characteristics [33], population dynamics [34], and biomass estimation [35,36]. However, the spatial distribution pattern of Picea schrenkiana in the TS is not clear, which has become the bottleneck for the scientific management and utilization of the plant resources. Therefore, we proposed an automatic extraction procedure for the spatial distribution of Picea schrenkiana based on GEE and the Jeffries–Matusita (JM) distance. The key research contents are as follows: (1) Can eliminating abnormal sample data improve the accuracy of land-cover classification in the TS of CA? (2) The difference in land-cover classification accuracy in the TS of CA under different seasonal images. (3) Feasibility analysis of using the JM distance to determine the importance of classification features. (4) Eight different scenarios were designed to determine the effects of different types of classification features and the best scenario was to extract the spatial distribution of Picea schrenkiana with a resolution of 10 m in the TS.

2. Materials and Methods

2.1. Study Area Description

This study was conducted in the TS (69°–95° E, 39°–46° N), with a total area of about 9.97 × 10⁵ km² (Figure 1). Due to its vast area and distance from the ocean, the TS is the farthest mountain system from the ocean in the world [37] and also the largest mountain system in the arid area of the world. The main topography of the TS consists of mountains, intermountain basins, and piedmont plains. The mean altitude of the ridge is about 4000 m and the highest peak has an altitude of about 7443.8 m [38]. Located in the northern temperate zone and far from the sea, the area is characterized by a temperate continental climate. The annual average temperature is about 5.0 °C [39], and the annual average precipitation ranges from 250 to 300 mm. The vegetation distribution of the TS has a typical vertical zone. From low altitude to high altitude, there are mountain steppe, meadow steppe, mountain coniferous forest, alpine steppe, alpine cushion plants, and snow-covered glaciers [33]. Picea schrenkiana is the dominant species in mountain conifer forests. According to the field survey data and existing literature (

Table 1. The main elevation range of Picea schrenkiana on the northern slope of the Tianshan Mountains.

Elevation	Reference
1500–2700 m	Wang, Ren et al. [34]
1600–2800 m	Li, Chang et al. [33]
1600–2800 m	Lan, Xiao et al. [40]
1500–2800 m	Luo, Xu et al. [35]
1500–2800 m	Li, Luo et al. [41]
1400–2800 m	Jiang, Zhu et al. [42]

), the elevation range of Picea schrenkiana on the northern slope of the TS in CA is from 1400 to 2800 m. To ensure an accurate extraction of Picea schrenkiana, a 1300 m contour was drawn from ArcGIS software and SRTM data with a spatial resolution of 90 m. The field enclosed by the contour at 1300 m is defined as the region studied in this paper.

2.2. Remote-Sensing Data

Sentinel-1 satellites are equipped with a dual-polarization C-band synthetic aperture radar (SAR) instrument, which provides SAR images with high temporal and spatial resolution. Sentinel-1 data are interferometric wide swath mode (IW) images, which is the standard mode over land [43] and corresponds to COPERNICUS/S1_GRD in GEE. The Sentinel-1 Toolbox has been used to preprocess the dataset, including thermal noise removal, radiometric calibration, and terrain correction [44]. Here, we exclusively utilized the vertical transmit/horizontal receive (VH) and vertical transmit/vertical receive (VV) bands. We also used Sentinel-2 data from the GEE archive. Sentinel-2 consists of two satellites (Sentinel-2A and Sentinel-2B), each carrying a multi-spectral imager (MSI) to acquire data in 13 spectral bands along a 290 km sun-synchronous orbit path every 5 days. Since Sentinel-2 MSI, Level-2A data were not provided fully on GEE’s data pool, Level-1C data were employed in this work, which corresponds to COPERNICUS/S2 in GEE. The datasets are products of the top of the atmosphere (TOA) after radiometric calibration and terrain correction.

Due to the wide coverage of the study area, it is necessary to combine images taken at different times to obtain high-quality image data covering the study area. At the same time, considering the factors of vegetation growth seasons, the images are divided into the four seasons of spring, summer, autumn, and winter (

Table 2. Sentinel-1/2 image selection schedule.

Seasons	Years	Months	Cloudy _Pixel_Percentage
Spring	2019–2020	4–6	30
Summer	2019–2020	7–9	20
Autumn	2019–2020	10–12	35
Winter	2019–2020	1–3	60

). The ee. Filter.calendarRange function was used to obtain images in the time range of the four seasons. The shapefile for the study area was imported into the GEE platform and the ee. filterBounds filter was employed to obtain images of the region of interest. For Sentinel-2 images, CLOUDY_PIXEL_PERCENTAGE of the four seasons was set to different values to remove all Sentinel-2 images of poor quality and obtain the number of pixel-by-pixel high-quality observation images in the study area (Table 2). The results showed that during the study period, the number of high-quality Sentinel-2 images from the four seasons in the study area ranged from 1 to 288 (Figure 2). The ee. Reducer.median function was then utilized to create a single image by calculating the median of all values on each pixel across the stack of all images. Finally, the Sentinel-1 composite image and high-quality Sentinel-2 composite image covering the study area were achieved for subsequent studies.

2.3. Training and Validation Sample Points

The accuracy of classification results is intimately related to the representativeness and quantity of samples [45]. Under realistic conditions for the TS and available land-cover data, the land-cover categories in this region were divided into seven types to separate Picea schrenkiana from the rest of the land uses. These include Picea schrenkiana, grassland, cropland, built, bare land, snow and ice, and waterbody. Due to the large extent of the study area, visual interpretation of a large number of sample points not only requires a lot of work, but also fails to properly distinguish between land-cover classes. To address this issue, we combined visual interpretation using Google Earth VHRI with semi-automatic extraction of sample points. On the one hand, the Picea schrenkiana samples were collected from VHRI on Google Earth. On the other hand, the non-Picea schrenkiana data (e.g., grassland, cropland, built, bare land, snow and ice, and waterbody) were collected by the semi-automatic extraction of sample points method. Specifically, semi-automatic extraction of sample points refers to the use of two tools (Create Random Points tool and Extract Multi Values to Points tool) in ArcGIS software and ESA World Cover data to obtain a large number of samples more accurately and quickly. Furthermore, the coordinates of 99 field survey Picea schrenkiana samples (Figure 1) were collected through GPS in the summer of 2019. Finally, a total of 8159 samples (

Table 3. The situation of Sample Collection.

Type	Picea schrenkiana	Grassland	Cropland	Built	Bare Land	Snow and Ice	Waterbody
Numbers	1524	3036	533	517	1704	275	570

) were collected, and 70% of these (5866) were used for model training and 30% (2293) for validation.

2.4. Classification Feature Input

Seven categories of features are applied for classification in this study (

Table 4. Summary of classification features applied in this research.

Categories	Features	Central Wavelength/Indices Formula	Data Source
Radar features	VV	vertical transmit/vertical receive	Sentinel-1
	VH	vertical transmit/horizontal receive	Sentinel-1
Spectral features	AEROS	443 nm	Sentinel-2
	BLUE	490 nm	Sentinel-2
	GREEN	560 nm	Sentinel-2
	RED	665 nm	Sentinel-2
	RDED1	705 nm	Sentinel-2
	RDED2	740 nm	Sentinel-2
	RDED3	783 nm	Sentinel-2
	NIR	842 nm	Sentinel-2
	RDED4	865 nm	Sentinel-2
	VAPOR	940 nm	Sentinel-2
	CIRRU	1375 nm	Sentinel-2
	SWIR1	1610 nm	Sentinel-2
	SWIR2	2190 nm	Sentinel-2
Red-edge features	NDVIre1	(RDED4 − RDED1)/(RDED4 + RDED1)	Sentinel-2
	NDVIre2	(RDED4 − RDED2)/(RDED4 + RDED2)	Sentinel-2
	NDVIre3	(RDED4 − RDED3)/(RDED4 + RDED3)	Sentinel-2
	NDre1	(RDED2 − RDED1)/(RDED2 + RDED1)	Sentinel-2
	NDre2	(RDED3 − RDED1)/(RDED3 + RDED1)	Sentinel-2
Spectral indices	NDVI	(NIR − RED)/(NIR + RED)	Sentinel-2
	EVI	2.5 × (NIR − RED)/(NIR + 6.0 × RED −7.5 × BLUE + 1)	Sentinel-2
	MNDVI	(RED − GREEN)/(RED + GREEN)	Sentinel-2
	NDWI	(GREEN − NIR)/(GREEN + NIR)	Sentinel-2
	LSWI	(NIR − SWIR1)/(NIR + SWIR1)	Sentinel-2
	MNDWI	(GREEN − SWIR1)/(GREEN + SWIR1)	Sentinel-2
	NDTI	(SWIR1 − SWIR2)/(SWIR1 + SWIR2)	Sentinel-2
	NDI45	(RDED1 − RED)/(RDED1 + RED)	Sentinel-2
Texture features	VV ASM	VV Angular Second Moment	Sentinel-1
	VV CON	VV Contrast	Sentinel-1
	VV CORR	VV Correlation	Sentinel-1
	VV SVAR	VV Sum Variance	Sentinel-1
	VV ENT	VV Entropy	Sentinel-1
	VH ASM	VH Angular Second Moment	Sentinel-1
	VHCON	VH Contrast	Sentinel-1
	VH CORR	VH Correlation	Sentinel-1
	VH SVAR	VH Sum Variance	Sentinel-1
	VH ENT	VH Entropy	Sentinel-1
Terrain features	ELE	ELEVATION	SRTM
	SLO	SLOPE	SRTM
	ASP	ASPECT	SRTM
	SHA	HILL SHADE	SRTM

). These include radar features, spectral features, red-edge features, spectral indices, texture features, and terrain features. Multiple lines of evidence showed that texture features, especially texture features based on SAR data, can avoid the phenomenon of “ same spectrum foreign objects” and “ same object foreign spectrum”, to improve the accuracy of remote-sensing image classification [46]. A gray-level co-occurrence matrix (GLCM), proposed by Haralick in 1973 [47], is the most common statistical analysis method for calculating texture features at present. Therefore, this study used the algorithm glcmTexture provided by GEE to compute the texture metrics with a window size of 3 × 3 and an offset of 1 × 1 from the GLCM around each pixel of Sentinel-1 bands (VV and VH) [48]. Considering that Picea schrenkiana is a single group tree species on the north slope of the TS at an altitude of 1400–2800 m, the ee. Algorithms. Terrain and ee. Image (“USGS/SRTMGL1_003”) function provided by the GEE platform was used to calculate terrain features.

2.5. The Method for Eliminating Abnormal Samples

The Box-plot developed by American statisticians in 1977 [49] is a common method to check whether there are outliers in a group of data. Compared with the Grubbs method, 3σ method, and Z-score method, the Box-plot has a wider application because it does not require the data to obey normal distribution to determine outliers. The 1st quartile (Q1), 2nd quartile (Q2), and 3rd quartile (Q3) are the numbers placed at 25%, 50%, and 75% of a set of data in ascending order, respectively [50] (Figure 3). The difference between Q1 and Q3 is called the interquartile range (IQR). The criteria for identifying outliers in the Box-plot are whether the data are greater than Q3 + 1.5IQR or less than Q1 − 1.5IQR in a group of data [51]. Minimum and maximum are the minimum and maximum values except for outliers. In the process of visual interpretation, Google Earth VHRI and semi-automatic extraction of sample points were used to determine land-cover categories. It is impossible to ensure that all sample points are pure and typical, which leads to the existence of outliers in the value of the extracted classification feature. For a specific sample, if the number of outliers in the 13 extraction values based on spectral features is greater than or equal to 2, the sample is considered to be an abnormal sample. Therefore, the outliers of the extracted values are identified based on the Box-plot and then the abnormal samples are judged and eliminated. By comparing the Kappa score and OA before and after eliminating abnormal samples, the influence of abnormal samples on classification accuracy was analyzed.

2.6. The Method for Classification Feature Optimization

Considering that some features selected in this study could not satisfy the requirement to distinguish Picea schrenkiana from other land-cover categories, redundant information existed in classification features [52], which affected the accuracy of mapping the spatial distribution of Picea schrenkiana. For this problem, the purpose of separability analysis is to evaluate the ability of each feature to distinguish Picea schrenkiana from other land-cover categories. There are three common methods to determine the separability between different categories: the separability index (SI) [53,54], the transformed divergence (TD) [55], and the JM distance [56,57]. Compared with the SI and TD, the JM distance can quantitatively evaluate the separability between different categories. Thus, the JM distance is used to evaluate the ability of each feature to distinguish Picea schrenkiana from other categories, to determine the best feature combination to extract Picea schrenkiana. The JM value is between 0 and 2, and its size represents the degree of separation of the seven land-cover categories under the selected features. When JM = 2, it indicates that seven types of land-cover categories are entirely separated under the selected classification features [58,59]. The expression for the JM distance is as follows:

$$\mathrm{JM}=2\left(1-{\mathrm{e}}^{-\mathrm{B}}\right)$$

(1)

where B represents the Bhattacharyya distance on a feature dimension. On the premise that normal distribution is satisfied, the Bhattacharyya distance between samples of different land-cover types is:

$$\mathrm{B}=\frac{1}{8}({\mathrm{e}}_1-{\mathrm{e}}_2{)}^2 \frac{2}{{\delta }_1^2+{\delta }_2^2}+\frac{1}{2}ln\left[\frac{{\delta }_1^2+{\delta }_2^2}{2{\delta }_1·{\delta }_2}\right]$$

(2)

where e_k represents the mean value of a certain type of feature; ${\delta }_{\mathrm{k}}^2$ represents the variance of a certain type of feature, where (k = 1, 2).

2.7. Random Forest-Based Classifier

RF algorithm is a non-parametric machine-learning algorithm of multiple decision-tree combinations proposed by Breiman Leo and Adele Cutler in 2001 [60]. In contrast to machine-learning algorithms such as decision trees and support vector machines, the RF algorithm is stable and effective in processing data multi-collinearity and dimensionality [61,62] without an over-fitting phenomenon [63]. Correspondingly, the RF algorithm was chosen and used in this study to map the spatial distribution of Picea schrenkiana in the TS of the CA.

The dominant construction procedure of the RF algorithm is as follows: based on the bootstrap sampling method, 2/3 of the samples are randomly chosen from the original sample set to construct the decision-tree model. Assuming that the features of each sample are M dimensions, m feature subsets (m << M) are randomly extracted from them, and the optimal feature is selected from these m features to split the node of the decision tree. The above procedure is iterated n times to obtain a random forest composed of n decision trees, and the optimal result is chosen from the n decision trees by voting for the classification result for each classification sample. In this study, a decision tree of size 50 was selected as optimal.

2.8. Description of Different Scenarios

Eight scenarios (

Table 5. Eight scenarios for mapping the spatial distribution of Picea schrenkiana.

Scenarios	Radar Features	Spectral Features	Red-Edge Features	Spectral Indices	Texture Features	Terrain Features	Features Optimization
S1		√
S2		√			√
S3		√		√
S4		√	√
S5	√	√
S6		√				√
S7							√
S8	√	√	√	√	√	√

) were designed to illustrate the effects of different types of classification features and the optimal classification features conducive to the extraction of Picea schrenkiana. We selected the best scenario and optimized the parameters of the classifier to extract the spatial distribution of Picea schrenkiana with a resolution of 10 m in the TS.

2.9. The Method for Accuracy Assessment

Accuracy evaluation is performed by constructing a confusion matrix based on the 2293 validation sample points and classification results from this study. The purpose of constructing the confusion matrix is to obtain accuracy assessment indexes, such as producer accuracy (PA), user accuracy (UA), the OA, and the Kappa score. The PA is the ratio of the number of pixels correctly assigned by the classifier to the actual total number of references in the category, which is used to measure the omission error of each category. The UA represents the ratio between the number of pixels correctly assigned to this category and the total number of pixels classified into this category by the classifier, which is used to measure the commission error of this category. The OA and Kappa scores are used to evaluate the overall classification accuracy of all categories. The four evaluation indices are calculated as follows:

$$\mathrm{PA}=\frac{N_{ij}}{N_{\ast j}}\times 100\%$$

(3)

$$\mathrm{UA}=\frac{N_{ij}}{N_{i\ast }}\times 100\%$$

(4)

$$\mathrm{OA}=\frac{\sum _{i=1}^q{N}_{ii}}{N}\times 100\%$$

(5)

$$\mathrm{Kappa}=\frac{[N\sum _{i=1}^q{N}_{ii}-{\sum }_{i=1}^q(N_{i\ast }\ast N_{\ast i})]}{[N^2-\sum _{i=1}^q(N_{i\ast }\ast N_{\ast i})]}$$

(6)

N_ij is the value of row i and column j in the confusion matrix, N_∗j is the sum of the j column of the confusion matrix, N_i∗ is the sum of row i of the confusion matrix, N is the total number of validation sample points, and q is the number of rows and columns of the confusion matrix.

3. Results

3.1. The Result of Eliminating Abnormal Samples

A total of 656 anomalous samples were selected by using the method described in Section 2.5. These include 117 Picea schrenkiana, 235 grassland, 56 cropland, 56 built, 139 bare lands, 24 snow and ice, and 29 waterbody. In this study, Picea schrenkiana was taken as an example to draw scatter plots of the extraction values of Picea schrenkiana samples in the spectral features (BLUE, GREEN, and RED) before and after the removal of abnormal samples. It can be seen from Figure 4 that after the removal of abnormal samples, the dispersion degree of the extraction values of Picea schrenkiana samples in the spectral features (BLUE, GREEN, and RED) decreased (Figure 4). In addition, the accuracy assessment index also changed. After the removal of abnormal samples, the OA and Kappa scores increased by 1.43% and 0.02, respectively, and the PA of Picea schrenkiana increased by 2.95%, but its UA decreased by 1.68% (

Table 6. Classification accuracy of abnormal samples before and after removal.

Type		Picea schrenkiana	Grassland	Cropland	Built	Bare Land	Snow and Ice	Waterbody	OA	Kappa
Before	PA/%	93.93	89.93	66.88	81.40	85.89	70.00	91.52	87.10%	0.83
	UA/%	95.43	84.52	79.85	90.32	83.47	78.87	96.79
After	PA/%	96.88	89.89	62.94	80.69	86.94	83.13	94.80	88.53%	0.85
	UA/%	93.75	87.65	81.08	81.82	85.78	90.79	95.35

Note: Darker red indicates lower values, while darker green indicates higher values.

).

3.2. Classification Results of Single Seasonal Images

The classification accuracy of different seasonal images is different; the OA and Kappa score was higher in spring and summer than in autumn and winter (

Table 7. Classification accuracy of single seasonal images.

Type		Picea schrenkiana	Grassland	Cropland	Built	Bare Land	Snow and Ice	Waterbody	OA	Kappa
Spring	PA/%	95.80	93.85	79.61	88.10	88.24	73.56	96.89	91.19%	0.88
	UA/%	96.71	89.91	90.98	93.28	87.88	81.01	97.50
Summer	PA/%	96.04	94.91	74.34	89.68	87.63	85.06	96.89	91.67%	0.89
	UA/%	96.26	88.52	91.87	91.13	91.53	90.24	98.73
Autumn	PA/%	93.94	91.48	75.00	92.06	84.79	70.11	87.58	88.36%	0.85
	UA/%	95.95	85.32	90.48	85.93	85.48	82.43	98.60
Winter	PA/%	95.34	93.14	69.74	87.30	78.09	62.07	80.75	86.39%	0.82
	UA/%	96.24	80.97	88.33	88.71	84.06	87.10	98.48

Note: Darker red indicates lower values, while darker green indicates higher values.

). The difference in the classification accuracy between different seasonal images is mainly caused by three land-cover categories: cropland, snow and ice, and waterbody. The main reason is that compared with other land-cover categories (Picea schrenkiana, grassland, built, and bare land), the seasonal variation in the three land-cover categories (cropland, snow and ice, and waterbody) is more obvious. The PA of Picea schrenkiana was the highest in summer and the lowest in winter (96.04% and 95.34, respectively), while the UA was the highest in spring and the lowest in winter (96.71% and 96.24%, respectively). Although the summer UA of Picea is lower than the spring UA, the summer PA of Picea schrenkiana increases by 0.24%–2.10% compared to the other three seasons. In addition, the OA and Kappa scores increased by 0.48%–5.28% and 0.01–0.07, respectively, in the summer compared to the other three seasons. As a result, the summer was used as the optimal sextant for mapping the spatial distribution of Picea schrenkiana. This conclusion is in accordance with the actual situation, as a great variety of plants are in a period of vigorous growth during the summer, and the physical properties and chlorophyll of the different plants are significantly different.

3.3. The Results of Feature Optimization

The JM distance of different classification features is considerably different (

Table 8. The JM distance of features and the accuracy of classification results of different numbers of features.

Feature	JM Distance	Kappa	Numbers
RED	1.01	0.41	1
SWIR2	0.94	0.608	2
GREEN	0.9	0.659	3
BLUE	0.89	0.73	4
RDED1	0.87	0.764	5
CIRRU	0.82	0.793	6
SWIR1	0.81	0.821	7
MNDVI	0.72	0.821	8
AEROS	0.7	0.828	9
VH ASM	0.68	0.829	10
VH ENT	0.6	0.836	11
LSWI	0.6	0.837	12
MNDWI	0.59	0.848	13
EVI	0.58	0.856	14
VAPOR	0.58	0.858	15
RDED2	0.56	0.864	16
NDTI	0.54	0.869	17
NIR	0.52	0.869	18
VV ASM	0.52	0.862	19
RDED4	0.51	0.868	20
RDED3	0.51	0.871	21
NDre2	0.49	0.873	22
NDre1	0.49	0.869	23
NDVIre1	0.47	0.87	24
NDVI	0.47	0.866	25
VV ENT	0.46	0.863	26
VH	0.39	0.873	27
NDVIre3	0.3	0.87	28
NDI45	0.3	0.869	29
SLOPE	0.28	0.886	30
VV	0.26	0.89	31
ELEVATION	0.24	0.894	32
VV CON	0.23	0.889	33
VV SVAR	0.22	0.89	34
NDVIre2	0.18	0.884	35
ASPECT	0.13	0.894	36
VH CON	0.11	0.892	37
VH SVAR	0.1	0.894	38
NDWI	0.09	0.89	39
HILL SHADE	0.07	0.894	40
VH CORR	0.02	0.89	41
VV CORR	0.02	0.898	42

Note: Darker red indicates lower values, while darker green indicates higher values.

). The feature with the largest JM distance is RED, while the smallest feature is VH CORR and VV CORR. The maximum distance is approximately 50 times the minimum distance. Based on the JM distances of the 42 categorical features in order from highest to lowest, the variation in the Kappa score was observed and the best combination of categorical features for mapping the spatial distribution of Picea schrenkiana was determined accordingly. With the increase in the number of features involved in classification, the Kappa score in the early stage (1–7 classification features) increased considerably, rapidly increasing from 0.410 to 0.821 (Table 8). In the middle period (8–17 classification features), the increased range of the Kappa score decreased significantly, showing a steady upward trend. In the later period (18–42 classification features), the Kappa score fluctuated slightly, fluctuating around 0.88. When the number of classification features used was 32 and 42, the Kappa score was 0.894 and 0.898, respectively. Considering that the difference between them is only 0.004, the first 32 classification features are selected as the final results of the feature optimization in this study.

The JM distances of the first 32 features were classified and analyzed (

Table 9. Jeffries–Matusita distance between different types of features.

Type of Feature	The Feature with the Highest Degree of Separation (Ranking)	The Feature with the Lowest Degree of Separation (Ranking)	The Average JM Distance
Spectral features	RED (1)	RDED3 (21)	0.74
Texture features	VH ASM (10)	VV ENT (26)	0.57
Spectral indices	MNDVI (8)	NDI45 (29)	0.54
Red-edge features	NDre2 (22)	NDVIre3 (28)	0.44
Radar features	VH (27)	VV (31)	0.32
Terrain features	SLOPE (30)	ELEVATION (32)	0.26

Note: Darker red indicates lower values, while darker green indicates higher values.

). The spectral features and texture features show high separation, while the radar features and terrain features show low separation. Among the spectral features, RED, GREEN, and BLUE have a high separation. The spectral features play an important role in forest mapping and are more prominent than other features in distinguishing Picea schrenkiana from other land-cover categories. In addition, the texture features play a key role in mapping the spatial distribution of Picea schrenkiana, as Picea schrenkiana interleaves with other land-cover categories such as grassland and bare land in patches.

3.4. Accuracy Analysis under Different Scenarios

Among the six scenarios (scenario 1- scenario 6), the PA and UA of mapping the spatial distribution of Picea schrenkiana were between 95.57%–97.44% and 95.40%–96.27%, respectively (

Table 10. Classification accuracy statistics under different scenarios.

Type		Picea schrenkiana	Grassland	Cropland	Built	Bare Land	Snow and Ice	Waterbody	OA	Kappa
Scenario 1	PA/%	96.27	92.90	61.84	75.4	86.00	81.61	95.03	88.75%	0.85
	UA/%	96.27	85.33	83.19	86.36	88.15	87.65	96.23
Scenario 2	PA/%	95.57	92.9	62.5	76.19	88.03	82.76	94.41	89.14%	0.86
	UA/%	95.79	86.36	85.59	84.96	87.50	91.14	96.82
Scenario 3	PA/%	96.50	94.08	65.79	75.40	87.63	81.61	95.03	89.84%	0.87
	UA/%	95.83	87.17	85.47	84.82	88.89	93.42	96.84
Scenario 4	PA/%	96.74	93.73	63.82	79.37	85.80	81.61	95.03	89.45%	0.86
	UA/%	95.62	86.46	88.18	82.64	89.62	87.65	96.23
Scenario 5	PA/%	96.74	92.54	66.45	78.57	87.02	80.46	95.03	89.36%	0.86
	UA/%	95.40	86.70	86.32	86.84	87.73	89.74	96.84
Scenario 6	PA/%	97.44	93.73	74.34	84.92	88.24	79.31	96.27	91.10%	0.88
	UA/%	95.87	89.19	90.40	89.17	88.41	94.52	97.48
Scenario 7	PA/%	96.74	95.50	73.68	91.27	87.63	83.91	95.65	91.93%	0.89
	UA/%	96.96	89.17	94.12	89.15	91.72	86.9	98.09
Scenario 8	PA/%	97.20	95.62	73.03	86.51	88.03	83.91	95.65	91.84%	0.89
	UA/%	96.53	88.79	93.28	88.62	91.75	91.25	98.72

Note: Darker red indicates lower values, while darker green indicates higher values.

). The OA of scenario 1 based on Sentinel-2 bands was 88.75%, and the Kappa score was 0.85. Based on the spectral features, the texture features (scenario 2), spectral indices (scenario 3), red-edge indices (scenario 4), radar features (scenario 5), and terrain features (scenario 6) were added, respectively, and the OA and Kappa scores were slightly increased. The OA scores increased by 0.39%, 1.09%, 0.7%, 0.61%, and 2.35%, while the Kappa scores increased by 0.01, 0.02, 0.01, 0.01, and 0.03, respectively. These results indicated that the spectral features had the best identification ability for Picea schrenkiana, while the other features had no obvious influence. However, the texture and terrain features still played a positive role in the land-use classification of the TS.

Compared with the other seven scenarios, scenario 7 obtained the best classification effect after optimizing all features by using the JM distance, and the OA and Kappa scores were 91.93% and 0.89 (Table 10), respectively. Compared to scenario 8, where all features participate in the classification, the OA and Kappa scores in scenario 7 are improved by 0.09% and 0.001, respectively, and the number of features in scenario 7 are less. Therefore, scenario 7 was selected as the final scheme for mapping the spatial distribution of Picea schrenkiana in the TS of CA.

3.5. Spatial Distribution of Picea schrenkiana in the Best Scenario

By comparing and analyzing different classification scenarios in Section 3.4, scenario 7 achieves the optimal classification effect, and its confusion matrix is shown in

Table 11. Confusion matrix of scenario 7.

	Picea schrenkiana	Grassland	Cropland	Built	Bare Land	Snow and Ice	Waterbody
Picea schrenkiana	415	13	0	0	0	0	0
Grassland	13	807	33	4	43	0	5
Cropland	0	5	112	2	0	0	0
Built	0	2	7	115	5	0	0
Bare Land	0	18	0	5	432	14	2
Snow and Ice	0	0	0	0	11	73	0
Waterbody	1	0	0	0	2	0	154
PA/%	96.74	95.5	73.68	91.27	87.63	83.91	95.65
UA/%	96.96	89.17	94.12	89.15	91.72	86.90	98.09
OA = 91.93%				Kappa = 0.89

Note: Darker red indicates lower values, while darker green indicates higher values.

. Picea schrenkiana and open water have the best mapping results, with both the PA and UA above 95%. The PA and UA of grassland, built, and bare land are both above 85%. The mapping accuracy of cropland and snow and ice is not high and the PA of cropland is only 73.68%. The spatial distribution of Picea schrenkiana in the optimal classification scenario (scenario 7) is shown in the figure and six places are randomly selected for detailed display. Picea schrenkiana is concentrated in the middle TS and scattered in the remaining areas (Figure 5). In addition, the pattern boundary of Picea schrenkiana based on the method presented in this paper has a strong consistency with the optical image (Figure 6).

A total of seven intervals with a separation of 300 m were set in the elevation distribution range of Picea schrenkiana, and the variation characteristics of Picea schrenkiana with altitude were observed. The distribution area of Picea schrenkiana first increases and then decreases with altitude (Figure 7, Elevation). Picea schrenkiana occupies the largest distribution area in the altitude range of 2100 to 2400 m, accounting for 29.7%. In total, 75.7% of the Picea schrenkiana is distributed in a range of elevation from 1800 to 2700 m. This study shows that Picea schrenkiana distributes in the TS of CA with a slope range of 0–69°. A total of seven slope ranges were set at intervals of 10° to observe the change in Picea schrenkiana distribution. In total, 33.9% of the Picea schrenkiana was distributed in the slope range of 30–40°, accounting for the largest proportion of all the intervals, and 66.4% of the Picea schrenkiana was distributed in the slope range of 20–40° (Figure 7, Slope). According to the classification method for aspects in forestry investigation, four aspects (semi-shady slope, shady slope, semi-sunny slope, and sunny slope) are set. Only 14.1% of the Picea schrenkiana were distributed in sunny and semi-sunny slopes, while 85.9% of the Picea schrenkiana were distributed in shady and semi-shady slopes (Figure 7, Aspect).

4. Discussion

4.1. The Influence of Feature Selection on Classification Results

In this study, the JM distance was used to explore the best classification features of land-use classification and map the spatial distribution of Picea schrenkiana in the TS. SWIR2, GREEN, BlUE, RDED1, CIRRU, SWIR1, MNDVI, AEROS, VH ASM, VH ENT, LSWI, MNDWI, EVI, and VAPOR were the top classification features in a sequence of importance (Table 8). Jiang et al. [42] selected an area of 10,000 km² in the TS for the remote-sensing identification of mountain coniferous forest and screened the ELEVATION, NDVI, spectral slope of red to near-infrared band, BLUE, RED, SWIR1. SLOPE is the most important classification feature of mountainous coniferous forest in the TS. SWIR1, RED, BLUE, SWIR1, and MNDVI (NDVI) in this study were consistent with Jiang’s research conclusion. The important reasons for BLUE and SWIR1 are that the bare land of the study area is large and the two classification features are sensitive to the soil background [64]. RED and NDVI are sensitive to vegetation and play a significant role in distinguishing vegetation from other land covers. ELEVATION and SLOPE, two categorical features, are of low importance in this study, which is inconsistent with the conclusions of Jiang et al. The main reason is that this study focuses on the importance of different features for mapping the spatial distribution of Picea schrenkiana. Therefore, only the JM distance between Picea schrenkiana and other land-cover types (grassland, cropland, built, bare land, snow and ice, and waterbody) was calculated. The distance between each of the two categories was not calculated. The immediate consequence of this practice would be to obscure the topographic features of Picea schrenkiana‘s distribution in the TS. In addition, the JM distance algorithm only evaluates the classification features based on the separation degree between different categories and does not consider the correlation between different features [65]. As a result, the classification features selected based on the JM distance selection may still have strong correlations, affecting the classification accuracy. At the same time, VH ASM and VH ENT are also of high importance in this study, because VH can obtain more vegetation structure information than VV [66].

4.2. The Spatial Distribution of Picea schrenkiana

The statistical results of this study showed that Picea schrenkiana was concentrated in the middle TS and scattered in other areas (Figure 5). Topographic factors such as elevation, slope, and aspect affect the spatial distribution of vegetation through the redistribution of vegetation habitat elements such as heat and water. The biomass of Picea schrenkiana showed a unimodal trend with the increase in altitude [35,67], and the biomass of Picea schrenkiana reached its peak at an altitude of 2100–2400 m, indicating that the altitude range was the most suitable for Picea schrenkiana to grow. The results of this study also show that the distribution area of Picea schrenkiana first increases and then decreases with the increase in altitude. The distribution area of Picea schrenkiana occupies the largest proportion in the altitude range of 2100 to 2400 m, accounting for 29.7% (Figure 7, Elevation). This is consistent with previous findings. Picea schrenkiana occupies the largest distribution area of 33.9% in the slope range of 30–40°, and 66.4% of the Picea schrenkiana is distributed in the slope range of 20–40° (Figure 7, Slope). In addition, 85.9% of the Picea schrenkiana is distributed in the shade slope and semi-shade slope (Figure 7, Aspect). These are consistent with the descriptions in the literature [68].

To sum up, although this paper has achieved good results in mapping the spatial distribution of Picea schrenkiana in the TS of CA (Figure 5), the impact of shadows was a major challenge for this work. In this area, the effect of removing shadows by terrain correction is not obvious, and the method for identifying the actual object type in the shadow according to the difference in the distribution position of shadows on the two images will have a large workload [42]. How to avoid the influence of shadows and choose the optimal feature selection algorithm [65] to more accurately identify the spatial distribution pattern of Picea schrenkiana in the TS is a problem worthy of attention in future studies.

5. Conclusions

The advantage of this work is being able to map the spatial distribution of Picea schrenkiana with 10 m resolution in the entire TS (9.97 × 10⁵ km²) by using GEE, the JM distance, the latest available satellite sensors, and a large number of training samples, overcoming the problem of insufficient computing capability for the large-scale monitoring of Picea schrenkiana. By using the Box-plot method to eliminate abnormal samples, the PA of Picea schrenkiana is increased by 2.95%, and the OA and Kappa scores of land-cover classification in the TS of CA are increased by 1.43% and 0.02. Compared to the other three seasons, the PA of Picea schrenkiana in summer increased by 0.24%–2.10%, the OA and Kappa scores of land-cover classification in the TS of CA by 0.48%–5.28% and 0.01–0.07. The separation obtained by the JM distance showed that spectral features and texture features played a key role in the mapping of Picea schrenkiana. The JM distance can seize the classification features that are most conducive to the mapping of Picea schrenkiana. The best classification features are RED, SWIR2, GREEN, BlUE, RDED1, CIRRU, SWIR1, MNDVI, AEROS, VH ASM, VH ENT, LSWI, MNDWI, EVI, VAPOR, RDED2, NDTI, NIR, VV ASM, RDED4, RDED3, NDre2, NDre1, NDVIre1, NDVI, VV ENT, VH, NDVIre3, NDI45, SLOPE, VV, and ELEVATION. Picea schrenkiana has a PA of 96.74%, and OA and Kappa scores of 91.93% and 0.89, respectively, for land-cover classification in the TS. Picea schrenkiana is concentrated in the middle TS and scattered in the remaining areas. In total, 85.7% of the Picea schrenkiana was distributed in the altitude range of 1500–2700 m, 66.4% of the Picea schrenkiana was distributed in the slope range of 20–40°, and 85.9% of the Picea schrenkiana was distributed in shady and semi-shady slopes.