An Automated Method for Surface Ice / Snow Mapping Based on Objects and Pixels from Landsat Imagery in a Mountainous Region

: Surface ice / snow is a vital resource and is sensitive to climate change in many parts of the world. The accurate and timely measurement of the spatial distribution of ice / snow is critical for managing water resources. Object-oriented and pixel-oriented methods often have some limitations due to the image segmentation scale, the determination of the optimal threshold and background heterogeneity. Therefore, this study proposes a method for automatically extracting large-scale surface ice / snow from Landsat series images, which takes advantage of the combination of image segmentation, the watershed algorithm and a series of ice / snow indices. We tested our novel method in three di ﬀ erent regions in the Karakoram Mountains, and the experimental results show that the produced ice / snow map obtained a user’s accuracy greater than 90%, a producer’s accuracy greater than 97%, an overall accuracy greater than 98% and a kappa coe ﬃ cient greater than 0.93. Comparing the extraction results under segmentation scales of 10, 15, 20 and 25, the user’s accuracy and producer’s accuracy from the proposed method are very similar, which indicates that the proposed method is more reliable and stable for extracting ice / snow objects than the object-oriented method. Due to the di ﬀ erent reﬂectivity values in the near-infrared band in the snow and water categories, the normalized di ﬀ erence forest snow index (NDFSI) is suitable for Landsat TM and ETM + images. This study can serve as a reliable, scientiﬁc reference for rapidly and accurately extracting ice / snow objects.


Introduction
Environmental changes and their impacts on natural ecosystems and human societies are topics of research in a wide range of scientific fields [1]. Surface ice/snow is among the most vital earth resources undergoing temporal and spatial changes as a consequence of climate change and other forms of environmental change in many parts of the world [2][3][4]. The economic, ecological and social effects of surface snow changes have been the subject of academic study for many years [5][6][7]. Furthermore, surface snow is one of the most important water resources of most rivers in high-mountain regions. It is not only the direct cause of spring flooding due to melting but also plays an important role in managing local water resources, preventing disasters, and forecasting and simulating snowmelt [8]. Due to global warming and corresponding glacier recession, ice/snow-covered areas in mountainous regions have decreased rapidly in recent decades. Therefore, monitoring and obtaining data on the dynamics of surface ice-snow in a timely manner are essential for the policy development process [9,10].
Since the Landsat Operational Land Imager (OLI) launched in 2013, many scholars have paid attention to temporal and spatial changes in the surface ice/snow [32]. However, few studies have focused on the differences among Landsat series images.
In summary, some problems in the previous methods, including segmentation parameters, classification thresholds and calculation speed, still exist [16,19,23,28,29,31]. Meanwhile, there is a lack of studies comparing the differences between ice/snow indices from Landsat TM, ETM+ and OLI images. Therefore, the objectives of this study are to (1) propose a method for integrating objects and pixels to lower the impacts of segmentation parameters and classification thresholds, and (2) to compare the adaptability of the ice/snow indices in Landsat series images.

Study Areas
The tested sites are located at the junction of China, Kyrgyzstan, Tajikistan, Afghanistan, Pakistan and India, spanning from 71 • E to 78 • E and from 34 • N to 40 • N (Figure 1), and the Karakoram Mountains pass through this region. The Karakoram Mountains have an average elevation of over 5500 m, a length of approximately 800 km and a width of approximately 240 km, possessing five mountain glaciers (Batura Glacier, Hispar Glacier, Biafo Glacier, Baltoro Glacier and Siachen Glacier) over 50 km and five peaks over 8 km (Chogori, Broad Peak, Gasherbrum I, Gasherbrum II and Teram Kangri) [33]. They are the most concentrated place for glaciers in the world outside the Arctic region. The study area is mainly semiarid and mostly belongs to a continental climate; the southern part of the study area belongs to a humid climate due to the influence of the Indian Ocean monsoon, but the northern part is extremely dry. The annual average precipitation does not exceed 100 mm at low altitudes, but it increases in snow-covered areas at high altitudes. At an altitude of approximately 5700 m, the average temperature in the warmest month is below 0 • C in this study area.
Remote Sens. 2020, 12,485 3 of 20 [31]. Since the Landsat Operational Land Imager (OLI) launched in 2013, many scholars have paid attention to temporal and spatial changes in the surface ice/snow [32]. However, few studies have focused on the differences among Landsat series images. In summary, some problems in the previous methods, including segmentation parameters, classification thresholds and calculation speed, still exist [16,19,23,28,29,31]. Meanwhile, there is a lack of studies comparing the differences between ice/snow indices from Landsat TM, ETM+ and OLI images. Therefore, the objectives of this study are to (1) propose a method for integrating objects and pixels to lower the impacts of segmentation parameters and classification thresholds, and (2) to compare the adaptability of the ice/snow indices in Landsat series images.

Study Areas
The tested sites are located at the junction of China, Kyrgyzstan, Tajikistan, Afghanistan, Pakistan and India, spanning from 71°E to 78°E and from 34°N to 40°N (Figure 1), and the Karakoram Mountains pass through this region. The Karakoram Mountains have an average elevation of over 5500 m, a length of approximately 800 km and a width of approximately 240 km, possessing five mountain glaciers (Batura Glacier, Hispar Glacier, Biafo Glacier, Baltoro Glacier and Siachen Glacier) over 50 km and five peaks over 8 km (Chogori, Broad Peak, Gasherbrum I, Gasherbrum II and Teram Kangri) [33]. They are the most concentrated place for glaciers in the world outside the Arctic region. The study area is mainly semiarid and mostly belongs to a continental climate; the southern part of the study area belongs to a humid climate due to the influence of the Indian Ocean monsoon, but the northern part is extremely dry. The annual average precipitation does not exceed 100 mm at low altitudes, but it increases in snow-covered areas at high altitudes. At an altitude of approximately 5700 m, the average temperature in the warmest month is below 0 °C in this study area.

Landsat Image
NASA has launched eight land satellites (Landsats) since 1972, and their products are widely used in various types of land and natural resource studies. Currently, only the Landsat 5, Landsat 7 and Landsat 8 products are available from the United States Geological Survey (USGS) Global Visualization Viewer (GLOVIS) portal. All the Landsat images are of product type L1T and have a scene quality score of nine, meaning perfect scenes with no detected errors. Furthermore, we selected  score of nine, meaning perfect scenes with no detected errors. Furthermore, we selected images from August to October to obtain perennial ice/snow. Detailed descriptions of the Landsat images are shown in Table 1.

Randolph Glacier Inventory
The Randolph Glacier Inventory (RGI) is a global inventory of glacier outlines, excluding the Greenland and Antarctic ice sheets, which was compiled at several institutions and can be downloaded from the website (http://www.glims.org/RGI/) [34]. Using simple and normalized band-ratio maps, most of the ice/snow outlines in the RGI were derived from various satellite imagery, including Landsat series images, ASTER IKONOS, and SPOT 5 [35]. RGI 1.0 was first released on 22 February 2012, and has been continuously updated since then. Currently, it has been updated to RGI 6.0, released on 28 July 2017. In this study, we selected RGI 4.0 (released on 1 December 2014) and RGI 5.0 (released on 20 July 2015) to indirectly verify the experimental results. According to their revision and acquisition times, RGI 4.0 matched Image C and Image F, and RGI 5.0 matched Image I.

Methods
Referring to some famous classification systems, including FAO LCCS and IGBP, we regard ice and snow as the same category [36]. The overall workflow has four steps ( Figure 2): (1) data preprocessing, which includes converting the digital number (DN) to top of atmosphere (TOA) reflectance data for every pixel and cloud mask, and calculating the NDSI, normalized difference forest snow index (NDFSI), normalized difference vegetation index (NDVI), NIR/SWIR and RED/SWIR; (2) applying a global threshold to extract ice/snow-covered areas after image segmentation using eCognition software (http://www.ecognition.com/suite), which is a famous remote-sensing classification tool; (3) marking ice/snow objects and building potential ice/snow-covered areas (unknown regions), and then using the watershed algorithm (https://docs.opencv.org/3.3.1/d3/db4/tutorial_py_watershed.html) to accurately delineate the boundary for each ice-snow object; and (4) evaluating the extraction results with the validation data and RGI datasets.

Data Preprocessing
The data preprocessing step includes three stages: converting the DN to the TOA reflectance, creating the cloud mask, and calculating the related indices.
The conversion of the raw pixel DN to the TOA reflectance is a basic step for Landsat image preprocessing. The DN value represents only the brightness characteristics of ground objects in the image, without any actual physical meaning. In contrast, the TOA reflectance is measured by spaceborne sensors flying outside the atmosphere, which can effectively reflect the spectral characteristics of ground objects and is widely used to calculate related indices. The conversion formula is given as follows from which we can extract gain and offset in the header files of the Landsat series images.

Data Preprocessing
The data preprocessing step includes three stages: converting the DN to the TOA reflectance, creating the cloud mask, and calculating the related indices.
The conversion of the raw pixel DN to the TOA reflectance is a basic step for Landsat image preprocessing. The DN value represents only the brightness characteristics of ground objects in the image, without any actual physical meaning. In contrast, the TOA reflectance is measured by spaceborne sensors flying outside the atmosphere, which can effectively reflect the spectral characteristics of ground objects and is widely used to calculate related indices. The conversion formula is given as follows from which we can extract gain and offset in the header files of the Landsat series images. Compared with other ground objects, clouds have high reflectivity in the SWIR band (2.08~2.35 µm). Based on this physical property of the Landsat images [37], we identified cloud-covered areas using single-threshold to reduce the influence of clouds on the experimental results. Furthermore, low cloudiness in the selected image also ensures the accuracy of experimental results.
Currently, the NDSI is widely used for large-scale surface ice/snow identification. The NDSI can effectively enhance the ice/snow information and suppress background information [9]. However, to some extent, using the NDSI to delineate snow-covered areas has some limitations due to heterogeneous false signals caused by the complexity of natural features. Some scholars have tried to construct a new ice/snow index according to the spectral characteristics of ice/snow in Landsat images, and introduced the NDVI and other indices to improve extraction results [27]. Although these studies can achieve good effects in their study areas, they are still limited to some extent by optical images from different sensors and different time phases. In this study, we identified snow-covered areas with different combinations of NDSI, NDFSI, NDVI, NIR/SWIR and RED/SWIR in the global extraction to obtain good extraction results (the spectral range of the SWIR band is 1.55-1.75 µm) [14,19,31]. The expressions of these indices are as follows RED/SWIR = P red P swir (6) where P green , P red , P nir and P swir represent the TOA values of the green band, red band, NIR band, and SWIR bands, respectively, for Landsat series images. Note that the spectral range for the SWIR band is 1.55-1.75 µm.

Global Extraction Based on Objects
Object-based classification, which is an effective method for extracting specific targets from Landsat series images, can avoid the noise of pixel-based classification [38]. Before object-based classification can be performed, image segmentation is a critical step. In this study, we adopted eCognition's multi-resolution segmentation algorithm, which is a bottom up region-merging technique starting with one-pixel objects. In numerous subsequent steps, smaller image objects are merged into bigger ones. Through this pair-wise cluster process, the weighted heterogeneity n h of the resulting image objects, where n is the size of a segment and h a parameter of heterogeneity, can be determined. In each step, that pair of adjacent image objects is merged, which results in the smallest growth exceeding the threshold defined by the scale parameter, causing the process to stop [39].
Heterogeneity in eCognition considers color and shape as primary objects. The increase in heterogeneity f has to be less than a certain threshold, which is also called a segmentation scale. f = w color * ∆h color + ∆h shape * w shape , w color ∈ [0, 1], w shape ∈ [0, 1], w color + w shape = 1 (7) ∆h color refers to difference in spectral heterogeneity, which allows multi-variant segmentation by adding a weight w c to the image channels c. with n merge number of pixels within merged object, n object_1 number of pixels in object 1, n object_1 number of pixels in object 2, σ c standard deviation within object of channel c. Subscripts merge refer to the merged object, and object 1 and object 2 prior to merge, respectively. ∆h shape refers to the difference in shape heterogeneity, which describes the improvement in the compactness and smoothness of an object's shape.
l is perimeter of object, b is perimeter of object's bounding box.
The DN values are multiples of the TOA values (the TOA values are from −1 to 1, and the DN values for Landsat TM are from 0 to 255), which makes σ c bigger, thereby increasing ∆h color and f, based on Formulas (7) and (8). Two different objects could only be merged when f is less than the segmentation scale. The increase in the heterogeneity of f could lower the probability of mixed-object and get better image segmentation results.
Therefore, the six bands (blue, green, red, nir, swir01 and swir02) from the original image were selected as the input layer to perform image segmentation using eCognition software, and then the related band ratios (NDSI, NDFSI, NDVI, NIR/SWIR, RED/SWIR) from the TOA image were applied to extract ice-snow information based on the segmentation results. We attempted to determine the segmentation scale several times, and the default values were used for w color and w shape (w compact and w smooth ) in this study. Due to the differences between different sensors in radiation resolution, the selected segmentation scales for Landsat TM, ETM+ and OLI were 15, 15, and 150, respectively.
A large ice/snow object can be divided into several small ones by image segmentation. These small ice/snow objects could be identified with a global threshold. In this study, we used two global thresholds to delineate snow-covered areas based on the object method ( Figure 2). Compared to the second global extraction, the threshold is larger in the first global extraction to ensure that the extracted ice/snow objects are accurate. The classification threshold (NDSI or NDFSI) is smaller in the second global extraction, to obtain potential snow-covered areas and provide support for building unknown regions.
Stratified random sampling was performed for different categories, and the calculations of the averages of the related indices for various categories in the study area are shown in Figure 3. According to Figure 3, there are obvious differences between the pure snow (i.e., non-shaded snow) and the other five classifications in terms of the NDFSI, NIR/SWIR and RED/SWIR for the Landsat TM and ETM+ images. However, the variances in NIR/SWIR and RED/SWIR are larger than those of the NDSI and NDFSI in pure snow sampling, as shown in Figure 3d. The NDSI has a good effect for extracting pure snow objects in Landsat OLI images. Therefore, we should extract pure snow mainly based on the NDFSI and NDSI, with NIR/SWIR and RED/SWIR as the auxiliary indices.
Ice/snow objects include two parts: pure snow and shaded snow. Compared to the pure snow type, shaded snow and water bodies have similar reflections in the five indices, which makes it difficult to directly extract shaded snow with a single criterion, as shown in Figure 3. We discovered that the differences between shaded snow and water bodies in the NDVI are larger than those in the other indices. Therefore, different criteria combinations for the five indices were applied to accurately extract ice/snow objects in this study.
We could set the classification threshold according to visual effects. Due to the differences in acquisition time and weather conditions for different images, the classification threshold varies from image to image. Taking image LE71480352002285SGS00 as an example, the criteria for the first global ice/snow extraction were as follows: (1)  more accurate than those in the second global extraction due to the strict criteria. However, there were some missing ice/snow objects in the first global extraction. Therefore, we lowered the criterion to include the potential ice/snow regions (also called unknown regions) and readjusted their boundaries in the local extraction based on pixels. other indices. Therefore, different criteria combinations for the five indices were applied to accurately extract ice/snow objects in this study.
We could set the classification threshold according to visual effects. Due to the differences in acquisition time and weather conditions for different images, the classification threshold varies from image to image. Taking image LE71480352002285SGS00 as an example, the criteria for the first global ice/snow extraction were as follows: (1) NDSI > 0.65, and (2) NDVI < −0.24 (Figure 4a). Based on the first global extraction results, we added a new criterion (NDSI > 0.58) for the second global ice/snow extraction (Figure 4b). The delineated boundaries of surface ice/snow in the first global extraction are more accurate than those in the second global extraction due to the strict criteria. However, there were some missing ice/snow objects in the first global extraction. Therefore, we lowered the criterion to include the potential ice/snow regions (also called unknown regions) and readjusted their boundaries in the local extraction based on pixels.

Local Extraction Based On Pixels
The object-based method regards the object as the basic classification unit, but the pixel-based method regards the pixel as the basic classification unit, which is their essential difference [24,25]. The delineated boundaries of surface ice/snow do not meet the mapping requirements, even after global extraction. Therefore, we needed to extract the unidentified ice/snow pixels from background information in the pixel layer.
As shown in Figure 4, the first depicted ice/snow areas are smaller than the actual ice/snow areas, while the secondary depicted ice/snow areas may be mixed with non-ice/non-snow pixels. To shrink

Local Extraction Based on Pixels
The object-based method regards the object as the basic classification unit, but the pixel-based method regards the pixel as the basic classification unit, which is their essential difference [24,25]. The delineated boundaries of surface ice/snow do not meet the mapping requirements, even after global extraction. Therefore, we needed to extract the unidentified ice/snow pixels from background information in the pixel layer.
As shown in Figure 4, the first depicted ice/snow areas are smaller than the actual ice/snow areas, while the secondary depicted ice/snow areas may be mixed with non-ice/non-snow pixels. To shrink the gaps between the automatically extracted ice/snow boundaries and the actual ice/snow boundaries, all the ice/snow bodies that were extracted in the first global extraction were numbered from 1 to n. Then, using the watershed algorithm, the boundary of every ice/snow body was readjusted within its potential ice/snow region.
The watershed algorithm regards the gradient image as a topographic surface and then floods from the minima based on region growing [40,41]. Finally, the image can be divided into watershed lines and catchment basins. Although the watershed algorithm has obvious advantages in terms of effectiveness and weak edge detection, over-segmentation cannot be avoided due to image noise. Therefore, a marker extraction method for the watershed algorithm was proposed [42,43]. The extracted ice/snow markers, representing the interior of different ice/snow bodies, are regarded as the minima of gray images and suppress all other gradient minima.
In the process of local extraction, we first numbered the depicted ice/snow bodies in the first global extraction, which were regarded as the sure foreground for the watershed algorithm (Figure 5a). The sure foreground, encircling a non-snow body, is segmented into multiple ice/snow bodies to ensure that the watershed algorithm works. Therefore, the segmented ice/snow objects need to be re-marked ( Figure 5b). Next, the buffer regions produced by the second global extraction results were spatially overlaid with the segmented ice/snow objects to develop the final markers (Figure 5c). The final markers include the sure ice/snow and unknown regions. Finally, based on the watershed algorithm, the boundary of every ice/snow object was automatedly modified within its unknown regions. The final extraction effect of ice/snow is shown in Figure 5d.

Accuracy Assessment
In this study, we adopted two methods for validating the experimental results from the qualitative and quantitative perspective. On the one hand, the extraction results and the RGI datasets were used to describe the boundaries of ice-snow from actual images and to show whether our method is effective from a qualitative perspective.
On the other hand, modifying the first global extraction results by visual interpretation were regarded as validation data ( Figure 2). We counted the corrected number of ice/snow pixels a , the commission number of ice/snow pixels b , the omission number of ice/snow pixels c and the corrected number of non-ice/snow pixels d . Six traditional indicators were selected-user's accuracy (UA), producer's accuracy (PA), commission error (CE), omission error (OE), overall accuracy (OA) and kappa coefficient from the confusion matrix [44], to evaluate the differences between the experimental results and the actual ice/snow boundaries for every image scene. The formulas of partial indicators are as follows

Accuracy Assessment
In this study, we adopted two methods for validating the experimental results from the qualitative and quantitative perspective. On the one hand, the extraction results and the RGI datasets were used to describe the boundaries of ice-snow from actual images and to show whether our method is effective from a qualitative perspective.
On the other hand, modifying the first global extraction results by visual interpretation were regarded as validation data (Figure 2). We counted the corrected number of ice/snow pixels a, the commission number of ice/snow pixels b, the omission number of ice/snow pixels c and the corrected number of non-ice/snow pixels d. Six traditional indicators were selected-user's accuracy (UA), producer's accuracy (PA), commission error (CE), omission error (OE), overall accuracy (OA) and kappa coefficient from the confusion matrix [44], to evaluate the differences between the experimental results and the actual ice/snow boundaries for every image scene. The formulas of partial indicators are as follows UA = a (a + b)

Surface Ice/Snow Areas Mapped Using Different Images
The developed method was implemented to identify the boundaries of surface ice/snow in the nine image scenes, as shown in Figure 6. We can see that the delineated boundaries of ice/snow from the developed method fit perfectly with the actual boundaries for different regions and time phases in this figure. Some weak boundaries were also effectively identified, and only a few ice-snow objects were missing. Although there were many shaded snow-covered areas in three images from region: 151035, they were still delineated with the proposed method.

Accuracy Assessment of The Extraction Results
To validate the extracted results for surface ice/snow, we modified the first global extraction results, converted them from vector to raster data, recorded the differences between the two results for each pixel and constructed a confusion matrix. A relatively straightforward comparison is shown in Table 2. Generally, large kappa coefficients (more than 0.93) and a large overall accuracy (greater than 98%) can be observed, which shows that there is great consistency between the experimental results and the reference results. Compared to the delineated boundaries from the proposed method, the depicted boundaries from the RGI datasets did not cover the snow area well. Although the snow-covered areas from the RGI datasets were reflected in the OLI images, there were some obvious omission errors and commission errors for snow. Moreover, the corrected areas from the RGI datasets were significantly less than the misclassified areas in Image F. From a qualitative perspective, the proposed method can meet the requirements of ice/snow mapping.

Accuracy Assessment of the Extraction Results
To validate the extracted results for surface ice/snow, we modified the first global extraction results, converted them from vector to raster data, recorded the differences between the two results for each pixel and constructed a confusion matrix. A relatively straightforward comparison is shown in Table 2. Generally, large kappa coefficients (more than 0.93) and a large overall accuracy (greater than 98%) can be observed, which shows that there is great consistency between the experimental results and the reference results.
For the ice/snow category, the user's accuracies ranged from 89.55% to 96.13%, and the mean accuracy was 93.70%; the producer's accuracies ranged from 98.01% to 99.26%, and the mean accuracy was 98.77%. These data intuitively show that the proposed method is suitable for extracting snow-covered areas. There is an interesting phenomenon that the user's accuracies are all lower than the producer's accuracies, which is mainly induced by some non-snow pixels that are misclassified as snow pixels. A global threshold was selected as a criterion for distinguishing ice/snow pixels from non-snow pixels during object-based global extraction. Although we used different combinations of the NDSI, NDFSI, NDVI, NIR/SWIR and RED/SWIR for different images to improve the global extraction results in this study, there were still some non-snow objects that were misidentified as ice/snow objects. The ice/snow pixels in unknown regions were extracted with the watershed algorithm to readjust the boundaries of the ice/snow objects during pixel-based local extraction, which reduced the omission error and improved the producer's accuracy for the ice/snow category. Moreover, it is demonstrated that the watershed algorithm has a good response to weak edges.
For the non-snow category, the user accuracies and producer accuracies are all higher than 97.9%, which is mainly because the non-snow-covered areas are larger than the ice/snow-covered areas in each image scene. Taking Image A as an example, the number of non-snow pixels (39,558,086) is greater than that of snow pixels (15,893,225) according to Table 3, which makes the user's and producer's accuracies of non-snow greater than 98%. Compared to the non-snow category, the user's and producer's accuracies of snow have obvious differences, due to the smaller covered areas, which is consistent with Zhao's research [45]. The number of omission pixels is only one-sixth of the commission pixels in Table 3, which indicates that the watershed algorithm can effectively identify ice/snow pixels from unknown regions in local extraction to improve the total accuracy.

The Influence of the Segmentation Scale on the Extracted Results
Image segmentation is a key step in the analysis of remote sensing imagery. The scale parameter is directly related to the quality and accuracy of object-oriented analysis in the segmentation process [26]. Selecting an inappropriate segmentation scale may cause mixed objects to merge or may fragment the ground objects, thereby decreasing the accuracy of the information extraction. Therefore, some researchers have quantified the difference between the actual segmentation results and the reference results with the potential segmentation error and Euclidean distance to evaluate the image segmentation results in a supervised manner [46,47]. They have also calculated the homogeneity within each ground object to evaluate it in an unsupervised manner [48]. However, the abovementioned methods are time consuming and labor intensive to implement, which is not suitable for large-scale mapping.
In this study, we determined the scale parameter based on visual effects to carry out image segmentation with eCognition in the global extraction phase. This method may not obtain an absolute optimal segmentation parameter, but it does not require sample-based image analysis and has high efficiency. Moreover, the boundaries of the ice/snow objects were readjusted with the watershed algorithm, which can eliminate the influence of the segmentation scale on the extraction results to some extent. As shown in Figure 7, Image H was segmented with different scales (10, 15, 20 and 25). The depicted boundaries under different segmentation scales are very similar for large ice/snow objects, and there are some gaps between the extraction results under different segmentation scales for small ice/snow objects. Some very small ice/snow objects are missed within a large segmentation scale because they will be mixed with other backgrounds. In other words, although there are different image segmentation effects under different scale parameters, the delineated boundaries for those are very similar.
According to Figure 8, the UA and PA of the proposed method (TPM) and the object-oriented method (TOM) change substantially when the segmentation scale varies from 25 to 200, illustrating that, regardless of the chosen method, the classification accuracy is affected by the segmentation scale. However, when the segmentation scale is less than or equal to 25, the user's accuracy, producer's accuracy, commission error, and omission error of TPM are obviously smaller than those of TOM. Under these scales (10, 15, 20 and 25), the maximum differences in the user's and producer's accuracies based on TPM are only 1.56 and 0.31, respectively; those based on TOM are 6.12 and 7.52, respectively. The results show that the extraction results from TOM are more dependent on imagery segmentation. Compared to TOM, TPM adds the local extraction based on the watershed algorithm, which readjusts the boundaries of ice and snow areas, thereby shrinking the gaps between the different extraction results within a certain segmentation scale.
For the ice/snow category, UA and PA of TPM are all higher than those of TOM, as shown in Figure 8. The average UA of TPM is 92.23% under these scales (10, 15, 20 and 25), but that of TOM is 87.90%. Similarly, the average PA of TPM is 97.94% under these scales (10, 15, 20 and 25), but that of TOM is 87.34%. The UA is similar for PA using TOM, and TPM has more obvious effects in improving PA than UA. This is mainly because TPM could reduce the omission number of ice/snow pixels c and increase the corrected number of ice/snow pixels a. Therefore, we can conclude that the proposed method can not only improve the extraction results, but also lower the impacts of segmentation scale to some extent.

The Influence of the Classification Threshold on the Extracted Results
The classification threshold is another important factor in remote sensing analysis. Using a relatively small threshold will increase the commission error, but using a relatively large threshold will increase the omission error. Therefore, selecting an optimal threshold is still a challenging task. Currently, researchers mainly extract the optimal threshold by sampling-based image analysis and expert experience [29,49]. However, it is difficult to distinguish ice/snow from the background effectively by using a simple threshold.
In this study, we used a relatively large NDSI or NDFSI value to extract pure ice/snow objects in the first global extraction, which could decrease the commission error but also increase the omission error. Then, based on the watershed algorithm, the total classification threshold for each ice/snow object was modified to readjust its boundary within its unknown region, which makes the classification threshold dynamic for each ice/snow object and decreases the omission error. As shown in Table 2, the omission errors for the ice/snow category are very small, ranging from 0.74% to 1.99%. Moreover, according to Figure 8, the omission errors of the proposed method are significantly fewer than those of the object-oriented method. These assessment results illustrate that the proposed method can effectively reduce the omission error and improve the overall accuracy.
There are also some similar "global-local" approaches for extracting specific ground objects [45]. The unknown regions are usually buffer zones directly produced by pure regions of a specific target, which can easily cause an increase in the commission error or omission error. As shown in Figure 4a,

The Influence of the Classification Threshold on the Extracted Results
The classification threshold is another important factor in remote sensing analysis. Using a relatively small threshold will increase the commission error, but using a relatively large threshold will increase the omission error. Therefore, selecting an optimal threshold is still a challenging task. Currently, researchers mainly extract the optimal threshold by sampling-based image analysis and expert experience [29,49]. However, it is difficult to distinguish ice/snow from the background effectively by using a simple threshold.
In this study, we used a relatively large NDSI or NDFSI value to extract pure ice/snow objects in the first global extraction, which could decrease the commission error but also increase the omission error. Then, based on the watershed algorithm, the total classification threshold for each ice/snow object was modified to readjust its boundary within its unknown region, which makes the classification threshold dynamic for each ice/snow object and decreases the omission error. As shown in Table 2, the omission errors for the ice/snow category are very small, ranging from 0.74% to 1.99%. Moreover, according to Figure 8, the omission errors of the proposed method are significantly fewer than those of the object-oriented method. These assessment results illustrate that the proposed method can effectively reduce the omission error and improve the overall accuracy.
There are also some similar "global-local" approaches for extracting specific ground objects [45]. The unknown regions are usually buffer zones directly produced by pure regions of a specific target, which can easily cause an increase in the commission error or omission error. As shown in Figure 4a, if the unknown regions were directly produced by ice/snow objects from the first global extraction, then the omission error would be obviously increased, thus reducing the overall accuracy. Compared to the traditional method, the unknown regions comprise two parts in the proposed method: (1) buffer zones produced by ice/snow objects from the first and second global extractions and (2) ice/snow objects from the second global extractions. The expanding method for unknown regions is more targeted, to reduce missed ice/snow objects and improve the overall accuracy.
Generally, the classification threshold has a great impact on the extraction results of the traditional method [29,50], which requires calculating the optical threshold to reduce the commission and omission errors. However, it is difficult to obtain the optimal classification threshold. A series of sampling-based image analyses may be necessary. Compared to the traditional method, there is no strict criterion for the optimal classification threshold in the proposed method, which is simple and effective. The classification threshold for each ice/snow object would be modified within its unknown region to meet its actual boundary. Therefore, the proposed method can decrease the impact of the classification threshold on the extraction results to some extent.

NDSI based on Landsat TM, ETM+ and OLI
Landsat images have been the most widely used data source for ice/snow mapping in recent decades, due to their long time series and high spatial resolution, although other images, including MODIS, Sentinel and SPOT5, have also been employed to map surface ice/snow [16,17,22,29]; the NDSI has also played an important role in this task. However, few studies have been conducted to compare the different responses of the NDSI to Landsat TM, ETM+ and OLI.
For Landsat TM and ETM+ images, ice/snow has a high reflectivity in the visible and NIR bands and has an absorption in the SWIR band [14,51]. The water body has similar spectral properties in the visible and SWIR bands, resulting in a similar reflectivity value for ice/snow and water bodies in the NDSI [27]. However, water bodies have a low reflectivity in the NIR band, which illustrates that the reflectivity of ice/snow is higher than that of water bodies in the NDFSI. As shown in Figure 3a,c, the NDSI makes no clear distinction between snow and water bodies. However, there is an obvious difference between them in terms of NDFSI. Therefore, the extraction results with the NDFSI are better than those with the NDSI in the Landsat TM and ETM+ images for the ice/snow category.
According to Figure 3b, the reflectivity of the NDSI in the Landsat OLI image has an obvious difference between the snow and water bodies, which is different from the Landsat TM and ETM+ images. Comparing these three sensors, we can discover that the spectral resolution from the Landsat OLI image is higher than that from the Landsat TM and ETM+ images, as shown in Table 4. The SWIR bandwidth (swir01) from Landsat OLI is only 0.1 µm, but those of Landsat TM and ETM+ are 0.2 µm. We can capture the subtle differences between different ground objects by improving the spectral resolution. In addition, the radiation resolution is increased from 8 bits to 16 bits, which enhances the contrast between different ground objects [9]. These changes produce obvious differences between surface ice/snow and water bodies in the NDSI from Landsat OLI images. Therefore, extracting ice/snow with the NDSI is suitable when using Landsat OLI images.

Conclusions
In the context of global warming, there is an urgent need to monitor earth surface ice/snow with Landsat series images in a timely and automatic manner. However, the accurate and rapid extraction of ice/snow in mountainous regions is still a challenging task due to many factors, such as image noise, image segmentation scale, classification criteria, mountain shadows and water bodies. Based on these challenges, a new method, that combines object-oriented classification and the watershed algorithm, was designed in this study to delineate the boundaries of ice/snow objects from Landsat TM, ETM+ and OLI images. The quantitative assessment of the extraction results by constructing a confusion matrix shows that the produced ice/snow map obtained a producer's accuracy greater than 97%, an overall accuracy greater than 98% and a kappa coefficient greater than 0.93. These results indicate that the developed method can obtain a good effect in extracting ice/snow objects in the mountainous region.
Comparing the extraction results under different segmentation scales and analyzing the impacts of the classification threshold, we discover that the proposed method is more stable and reliable for delineating the boundaries of ice/snow compared to the object-oriented method. This method, based on objects and pixels, can be used not only to extract ice/snow, but also to extract water bodies and other single objects, and it has broad applicability. Furthermore, due to the differences in the spectral and radiation resolutions, different indices and their combinations should be applied during global extraction when using Landsat series images.
Author Contributions: X.W. drafted the manuscript and was responsible for data preparation, experiment, and analysis. X.G. participated in data collection. X.Z. was responsible for the research design and analysis, designed and reviewed the manuscript. W.W. and F.Y. validated the experimental results. All authors contributed to the editing and reviewing of the manuscript. All authors have read and agreed to the published version of the manuscript.