IMPROVED ON SNOW COVER EXTRACTION IN MOUNTAINOUS AREAS BASED ON MULTI-FACTOR NDSI DYNAMIC THRESHOLD

Snow cover is one of the most active elements of the cryosphere and plays an important role in the surface radiation budget and water balance. Optical satellite remote sensing has become an important tool for snow identification and monitoring. The Sentinel-2 A/B satellite has become an important data source for snow cover extraction because of its high radiation resolution, which reduces the problem of snow/ice saturation in remote sensing observations. The normalized differential snow cover index (NDSI) and the snow cover index considering the forest cover area (NDFSI) are important methods for snow cover extraction, but due to the strong spatial heterogeneity and fast change speed of snow cover in mountainous areas, using the classical fixed threshold method to extract snow cover will result in a large omission error. In this paper, a dynamic threshold method is constructed by synthesizing the effects of the type of snow cover period, aspect and snow underlying land coverage type on the snow cover NDSI/NDFSI. Compared with the high-resolution GF-2 snow map results, the dynamic threshold method has higher accuracy in extracting snow cover, and the overall classification accuracy, omission error, commission error and Kappa coefficient are 97.70%, 0.20%, 11.17% and 0.93 respectively. The dynamic threshold method is used to extract snow cover in a snow cover period in the Babao River Basin. The snow cover rate in the basin fluctuates greatly with time, and the spatial distribution of snow cover is uneven, with more snow cover in mountainous areas and rapid changes in snow cover in river valleys. * Corresponding author E-mail address: zyl0322@nwnu.edu.cn (Y. Zhang).


INTRODUCTION
Snow is an important type of land cover, which plays an important role in surface radiation budget and water balance, and has an important impact on the natural environment and human social activities (Rosenthal and Dozier, 1996;Cohen and Entekhabi, 2001;Henderson et al., 2018). On the other hand, snow meltwater is an important fresh water for ecosystems in arid and semi-arid regions (Pulliainen 2006, Li et al., 2019, and snow cover area (SCA) is a key input parameter for hydrological cycle and climate studies. Therefore, how to accurately extract snow cover information has undoubtedly become an important research content in earth system science.
Traditional snow monitoring methods at meteorological stations are time-consuming and labor-intensive, while satellite remote sensing with large-area synchronous observations has become an important tool for snow monitoring (Negi et al., 2009). Optical satellite remote sensing has been widely used in snow information extraction research. Because snow cover has high reflectance in the visible light band and low reflectance in the short-wave infrared band, whereas clouds have high reflectance in the above two bands. Dozier (1989) proposed the concept of normalized difference snow index (NDSI), which can effectively distinguish snow from clouds and other ground objects. Hall et al., (1995) developed the SNOMAP algorithm for snow recognition based on NDSI, which has been successfully applied to the global snow cover products of the American Moderate Resolution Imaging Spectroradiometer (MODIS). Wang (1999) found that the SNOMAP algorithm based on NDSI is still the best snow cover extraction method for optical remote sensing. In addition, many scholars have found that NDSI cannot extract snow cover well when the snow cover underlying land type is woodland. Wang et al., (2015) conducted spectral analysis on forest regions with and without snow-covered, and proposed a normalized difference forest snow index (NDFSI), the threshold is set to 0.4 to extract the snow-covered forest, which improves the snow cover recognition accuracy.
The NDSI threshold is the key to snow cover extraction. MODIS Global Snow Products released by the National Ice Snow Data Centre (NSIDC) set NDSI threshold as 0.4 (Hall, 1998), which was later adopted by many domestic and foreign scholars. However, many researchers found that the 0.4 threshold was not suitable for all regions, and the optimal threshold of NDSI was affected by the satellite sensor types, local terrain conditions, snow characteristics and underlying land types (Härer et al., 2018;Zhang et al., 2019;Hao, 2021). The distribution of snow cover in mountainous areas is complex, with strong spatial heterogeneity and variation. Therefore, the use of the classical fixed threshold method to extract snow cover will result in a large commission error. The Sentinel-2 A/B satellite has become an important data source for snow extraction because of its high radiation resolution, which reduces the problem of snow/ice saturation in remote sensing observations (Kääb et al., 2016;Zhang et al., 2020).
Using a fixed threshold to extract mountain snow cover has low accuracy. This paper proposes a dynamic NDSI/NDFSI threshold mountain snow cover extraction method based on Sentinel-2 images, so as improve the snow cover extraction accuracy in mountain areas.

Study area overview
The Heihe River Basin is the second largest inland river basin in China, the Babao River Basin located in the southeast of the Heihe River Basin (37°43′15″N-38°18′24″N, 100°00′35″E-101°09′04″E), as shown in Figure 1. The terrain in the basin is complex, with altitudes ranging from 2603 m to 4969 m, with an average altitude of 3604 m. The Babao River Basin is surrounded by high mountains, and snowfalls frequently in winter, with more snowfall from October to April of the following year. The annual average temperature in the basin is -1°C, and the glaciers and snow cover areas are widely distributed, and there are glaciers and snow belts above 4500 m (Ge et al., 2015). The multi-year snow cover rate is 0.48%, and the seasonal snow cover rate is 17.53%, which belongs to stable seasonal alpine snow cover .

Figure 1.
Overview of the study area.

Data source and preprocessing
There are five main data sources used in the study, as shown in Table 1. 1) Sentinel-2 L1C products are ortho-corrected top of the atmosphere (TOA) spectral reflectance products. In this paper, a total of 93 scene images of 32 phases in a snow cover period in the study area were screened for simultaneous atmospheric and topographic correction, that is, topographic normalization; 2) ALOS 12.5 m digital elevation model data (DEM) is used for auxiliary data of L1C product topographic normalization processing and generating terrain factors; 3) GF-2 data is used to verify the accuracy of snow cover extraction; 4) Google Earth images are used for auxiliary verification data of snow cover extraction accuracy in forest areas; 5) The surface cover dataset ESA_WorldCover is used to extract the type of snow underlying land cover in the study area.  Figure 2, the snow cover extraction in mountainous areas includes three steps:1) First, this study uses the Sen2Cor topographic normalization tool provided by ESA to obtain a user-defined Sen2Cor-L2A (S2C-L2A) data by configuring auxiliary data and activating topographic correction parameters to remove the influence of terrain shadows on snow cover extraction. 2) The NDSI/NDFSI dynamic threshold method was constructed by integrating the snow cover period, aspect and snow cover underlying land type. 3) Finally, the dynamic NDSI/NDFSI threshold are used to extract snow cover in rugged terrain, and analyze the spatial and temporal distribution characteristics of snow cover during a snow cover period in the Babao River Basin.  Hall (1995) developed the SNOMAP algorithm for snow cover extraction, and NDSI is the core content of the SNOMAP algorithm. In the California TM data test, when 50% of the pixels are covered by snow cover, the NDSI value of the corresponding TM image pixel is greater than or equal to 0.4 (Hall et al., 2002). When using NDSI to extract snow cover, it is found that the NDSI value of some water bodies is also greater than or equal to 0.4. Therefore, the SNOMAP algorithm adds a discriminant condition, and the reflectance in the near-infrared band is set higher than 0.11 to exclude water bodies. In addition, many dark matter also have high NDSI values, and the SNOMAP algorithm also supplements the green band greater than or equal to 0.1 as a discriminant condition to remove the influence of dark matter on snow cover extraction. Finally, the pixels satisfying the above three conditions at the same time are determined as snow cover, and the complete expression of the SNOMAP algorithm is shown in formula (1). Due to the computational simplicity and effectiveness of the NDSI index in extracting snow cover, this study retains the use of NDSI, but its threshold size is the key to snow cover extraction. The spatial heterogeneity of snow cover in mountainous areas is strong. The SNOMAP algorithm uses a fixed NDSI threshold to extract snow cover, which ignores the temporal changes of snow spectral information. The NDSI threshold of mountain snow cover should be constructed with dynamic threshold to improve the accuracy of mountain snow cover extraction.

Identification method of snow cover in non-forest area and forest area
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France In the formula, 3  , 8a  , and 11  are the reflectance of the 3rd, 8ath, and 11 bands of Sentinel-2, respectively.
For the snow cover in the forest area, this study uses the NDFSI index to extract, and the NDFSI index is calculated as shown in formula (2). Compared with other methods, the NDFSI index method is simpler and more effective. However, when the existing research uses a fixed NDFSI threshold for forest area snow cover extraction, the effects of forest coverage, forest canopy structure and density on the spectral reflectance of forest area snow cover are ignored.
8a 11 In the formula, 8a  , and 11  are the reflectance of the 8a band and the 11th band of Sentinel-2, respectively.

Dynamic threshold construction method
Before constructing NDSI/NDFSI dynamic threshold, it is necessary to analyze the influence of each factor on snow NDSI/NDFSI. According to previous research results, the influence of snow cover period, snow underlying land coverage type and aspect on NDSI/NDFSI was analyzed. When analyzing the effects of the three factors on snow NDSI/NDFSI, a large number of NDSI/NDFSI samples need to be randomly collected for different aspects, snow underlying land coverage types and different snow cover period stages. In order to determine that the selected samples are all NDSI/NDFSI corresponding to snow cover, this study uses the Sentinel-2 false color image and Google Earth images as a reference to select samples. In order to increase the representativeness of the samples, 2000 samples were selected on the images with no cloud and less cloud cover in the accumulation period, the stable period and the ablation period, respectively.
Based on the analysis of snow cover period, snow underlying land coverage type and aspect, the different effects of the three factors on snow cover NDSI/NDFSI are divided into different categories. The snow cover period is divided into accumulation period (AC), stable period (ST) and ablation period (AB), snow cover underlying land type are divided into grassland, woodland and other three types, and the aspect is divided into three types: shady aspect/semi-shady aspect, sunny aspect/semi-sunny aspect and flat. Based on the two factors of snow cover underlying land type and aspect, the threshold are divided into seven types: grassland sunny aspect (GSU), grassland shady aspect (GSH), woodland sunny aspect (WSU), woodland shady aspect (WSH), other sunny aspect (OSU), other shady aspect (OSH) and FLAT. In order to make the obtained threshold reliable, through image screening, images of three phases with better snow cover distribution are selected as threshold determination images in three stages of snow cover period, and the final threshold is the average value of the three phases.
Since the spatial resolution of the Sentinel-2 satellite is 10 m, which can reflect the snow situation on the ground, this paper uses the visually interpreted snow pixels as the "true value". When determining the optimal threshold, the Sentinel-2 snow map under different NDSI/NDFSI threshold conditions is obtained by continuously changing the NDSI/NDFSI threshold and increasing it with a step size of 0.01. When the ratio of Sentinel-2 snow map to "true" snow pixel is closest to "1", the threshold at this time is considered to be the optimal threshold.

Snow cover extraction accuracy evaluation method
Since the spatial resolution of Sentinel-2 images is 10 m, and the spatial resolution of GF-2 images after fusion is 0.8 m, according to the snow cover verification method proposed by Hall et al. (1995), this paper will use the MeanVis method to extract GF-2 snow map is used as the "true value" to verify the accuracy of the snow cover extracted by the dynamic threshold method. Based on the confusion matrix, the overall classification accuracy (OA), commission error (CE), omission error (OE) and Kappa coefficient were used to quantitatively evaluate the snow cover accuracy extracted using dynamic threshold. The overall classification accuracy is the ratio of the number of correctly classified samples to the total number of samples. The commission error is the ratio of the pixels classified as snow but actually non-snow to the non-snow pixels in the classification result. The omission error is the ratio of the pixels classified as snow-free but the actual snow to the snow pixels in the classification result. The Kappa coefficient is used for consistency check and is also an indicator of classification accuracy.  Table 2. Definition of confusion matrix and evaluation index.

Dynamic threshold construction results
By constantly changing the thresholds of NDSI/NDFSI and comparing with the "true value" snow map, the threshold closest to the "true value" snow map is selected as the optimal threshold. Determine the optimal NDSI threshold for snow cover in non-forest areas, and determine the optimal NDFSI for snow cover in forest areas. Figure 3 shows the process of determining dynamic threshold during the accumulation period on November 06, 2020, where the larger circle represents the optimal thresholds. It can be seen from the figure that with the increase of the threshold, the fewer snow pixels are extracted, and the ratio of the snow pixels to the "true value" snow map gradually decreases. And it can be found that the NDSI thresholds of shady aspect of different snow underlying land types are larger than those of sunny aspect , which is mainly caused by the distribution of snow on shady aspect. The magnitudes of the seven threshold types are also quite different. The final optimal thresholds for GSU, GSH, WSU, WSH, OSU, OSH and FLAT were 0.14, 0.29, 0.31, 0.38, 0.21, 0.31 and 0.35, respectively.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France Figure 3. Distribution of snow pixel ratio with thresholds in sentinel-2 snow map and "true" snow map.
According to the above process, the NDSI/NDFSI thresholds of the other eight phases were determined respectively. The optimal NDSI/NDFSI thresholds for the final accumulation period, stable period and ablation period were the average values of the three phases, as shown in Table 3. It can be seen from the table that the NDSI/NDFSI threshold of the seven types are different, and the threshold in the accumulation period and the ablation period are both smaller than those in the stable period. The NDSI thresholds of snow-covered non-forest in the three stages of the snow cover period are all smaller than the standard fixed threshold of 0.4, which is consistent with the research results of Hao et al. (2008). The small threshold is mainly due to the rapid ablation and change of snow cover in the Qinghai-Tibet Plateau, and the small snow thickness leads to a generally small threshold. Compared with the NDFSI threshold of snow-covered forest of 0.4 proposed by Wang et al. (2015), this paper determined that the NDFSI thresholds of different stages of snow cover are less than 0.4. This also shows that using a fixed threshold of 0.4 will result in a large amount of snow cover that cannot be extracted.  Table 3. Dynamic threshold of snow cover period.

Dynamic threshold extraction snow cover accuracy verification
On March 23, 2020, Sentinel-2 and GF-2 transited part of the Babao River Basin at the same time. The snow cover of the GF-2 image extracted by the MeanVis method is used as the "true value" snow map as shown in Figure 4(b). Using the dynamic threshold method to extract seven types of snow cover is shown in Figure 4(d). Compared with the false color images, the dynamic threshold method can more accurately extract seven types of snow cover. And it can be seen that the snow cover on shady slopes is always more than that on sunny aspect in different snow underlying land coverage types, and the snow cover is mostly distributed on grass and flat ground. Figure 4(e) shows the snow cover results of the whole image extracted by the dynamic threshold method. Compared with the GF-2 snow map, the snow cover extracted by the dynamic threshold method is similar, and the visual results are better.   Table 4. The accuracy of snow cover by the dynamic threshold method.

Analysis of dynamic change characteristics of snow cover area
Through image screening, Sentinel-2 images of 32 phases of a snow cover period from 2020 to 2021 were selected, of which the accumulation period was 13 phases, the stable period was 8 phases, and the ablation period was 11 phases. The snow cover in a snow cover period is extracted by the dynamic threshold method, and its temporal variation and spatial distribution characteristics are analyzed respectively. Figure 6 shows the spatial distribution of snow cover extracted by the dynamic threshold method. It can be seen from the figure that the spatial distribution of snow cover during the accumulation period in the Babao River Basin increases with time, the stable snow cover is mostly scattered in the mountainous areas with higher altitudes, and the snow cover in the valley area varies greatly. It can also be seen from the figure that the snow cover in October is less distributed in the accumulation period, and is only scattered in high-altitude mountainous areas. Since November, the snow cover has continued to increase. Large-scale snowfall occurred after November 26, and the snow cover remained stable in a wide range in December. There is large-scale snow cover in highaltitude mountains and river valleys. Figure 7 shows the spatial distribution of snow cover in 8 phases during the stable period. It can be seen from the figure that except for the large-scale snowfall on January 25, 2021 in the study area, the distribution of snow cover in other periods is relatively small. Different from the accumulation period, the snow cover in the stable period is sporadically distributed outside the mountains with higher altitudes, and the river valleys are also sporadically distributed. Compared with the accumulation period, the distribution of snow cover in the stable period is relatively small. Figure 8 shows the spatial distribution of snow cover in 11 phases during the ablation period . It can be seen from the figure that the spatial distribution of snow cover varies greatly during the ablation period. Moreover, it can be seen that the snow cover during the ablation period is not clear with the terrain texture due to the rapid ablation due to the increase of temperature. During the whole ablation period, except for the large-scale snow cover distribution on March 1 and April 15, the valley area was not covered with snow cover during the rest of the ablation period, and the snow cover was only distributed in high-altitude mountainous areas. In general, the spatial distribution of snow cover in the Babao River Basin varies greatly, and the spatial distribution of snow cover at different stages of the snow cover period is different. Figure 5 shows the variation of snow cover ratio with time in the Babao River Basin. It can be seen from the figure that the snow cover rate of Babao River fluctuates greatly with time, the maximum snow cover rate reaches 86.58%, and the minimum snow cover rate is only 4.66%, the difference between the two is 81.92%. And the snow cover rate during the accumulation period did not continue to increase with time, but showed a trend of first increasing and then decreasing, and the overall snow cover rate was increasing; In the stable period, the snow cover rate did not stabilize at the maximum snow cover rate, but only stabilized at 15% to 20%; Similar to the snow cover rate in the accumulation period, the fluctuation in the ablation period was larger, and there was no continuous decrease, and the change trend was first decreased, then increased and then decreased. The above situations are all related to the fast melting and changing characteristics of snow cover on the Qinghai-Tibet Plateau.

Error analysis of dynamic threshold snow cover extraction
Figure 9(b) shows the error map of snow cover extraction by the dynamic threshold method. It can be seen from the figure that the error of snow cover extraction by the dynamic threshold method is mainly concentrated in the transition zone between the snow-covered area and the snow-free area. This part of the snow cover is also a difficult point for the existing extraction snow cover algorithm. Because most of the pixels in these locations are mixed pixels, their properties are unstable, and they are more easily affected by external environmental factors such as light conditions and atmospheric conditions, resulting in drastic changes in NDSI, which increases the error of snow cover extraction. Although the dynamic threshold method is expected to improve the snow cover extraction accuracy compared with the fixed threshold , it can be seen from the figure that the dynamic threshold does not extract the snow cover in the transition zone between the snow-covered area and the snow-free area. This is also the direction of future improvement.

Comparison of snow cover extraction accuracy between dynamic threshold and fixed threshold
Taking the image of April 5, 2021 as an example, the spatial differences of the snow cover extraction results with dynamic threshold and fixed threshold were obtained ( Figure 10). From the difference map of snow cover extraction results, it can be seen that the differences between dynamic threshold and fixed threshold extraction snow cover are mainly concentrated in three areas. The first is the snow-covered forest area as shown in Figure 10(a); the second is the transition zone between the snow-covered and snow-free areas as shown in Figure 10(b); the third is the thin snow cover area as shown in Figure 10(c). In the experiment, it was found that the difference caused by the transition zone between the snow-covered area and the snowfree area is due to the thin snow cover thickness. The pixels in these positions are mostly mixed pixels, and the snow distribution is seriously trivial. Extracting this part of the snow cover more difficult. Therefore, the transition area and thin snow area are classified into one category for difference analysis. The forest area and thin snow cover area with large differences in snow cover extracted by dynamic threshold and fixed threshold are partially enlarged and analyzed. Figure 11 shows the results of dynamic threshold and fixed threshold extraction snow cover in the thin snow cover area. It can be seen from the rectangular frame in the figure that in the area with relatively thin snow cover, dynamic threshold can identify a large amount snow cover while fixed threshold cannot, and more fragmented snow cover is extracted, indicating that the use of dynamic threshold can eliminate the extraction error caused by a single threshold to a certain extent, so that snow cover can be extracted more accurately. Figure 11. Snow cover extraction results in thin snow cover areas: (a) Sentinel-2 false color image; (b) dynamic threshold extraction snow cover results; (c) fixed threshold extraction snow cover results Figure 12 shows the results of extracting snow cover with dynamic threshold and fixed thresholds in the forest area. It can be seen from the rectangular frame in the figure that dynamic threshold can accurately extract the snow cover in the forest area while fixed threshold cannot, and the phenomenon of snow cover leakage is serious. This is mainly due to the fact that when the forest coverage is high and the snow cover is almost entirely below the canopy, the vegetation canopy will greatly weaken the spectral signal reflected by the surface snow cover, resulting in a decrease in the NDFSI threshold. However, when the NDFSI threshold is fixed for snow cover extraction in forest areas, the selected threshold is too high, so that snow cover in forest areas cannot be extracted. The NDFSI dynamic threshold method reduces the NDFSI threshold according to the forest coverage and aspect, which overcomes the insufficiency of the fixed threshold.

CONCLUSION
The distribution of snow cover in mountainous areas is complex, with strong spatial heterogeneity and fast change speed. Using the classic fixed threshold method to extract snow cover will result in a large omission error. Therefore, in order to improve the extraction accuracy of snow cover in mountainous areas, NDSI/NDFSI dynamic threshold were constructed by integrating the type of snow cover period, aspect and snow underlying land coverage type. At the same time, the dynamic threshold determined in this paper is used to extract the snow cover in a snow cover period (2020.10-2021.05) in the Babao River Basin, and analyze its dynamic changes in time and space.
(a) (c) The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France The NDSI/NDFSI thresholds were divided into seven types by synthesizing the effects of snow cover period, aspect and snow underlying land coverage type on snow cover NDSI/NDFSI. The NDSI/NDFSI thresholds of GSU, GSH, WSU, WSH, OSU, OSH and FLAT in three stages of snow cover period were determined respectively. Compared with the standard NDSI/NDFSI threshold, the dynamic threshold determined in this study is smaller, which means that snow pixels cannot be identified by using a fixed threshold to extract snow cover. Compared with the GF-2 snow map, the OA, OE, CE and Kappa coefficients of the snow cover extracted by the dynamic threshold method are 97.70%, 0.20%, 11.17% and 0.93, respectively, and the snow cover extraction accuracy is higher. The dynamic threshold method proposed in this paper can be applied to the snow cover extraction in other mountainous areas.
The snow cover rate in the Babao River Basin fluctuates greatly with time. The snow cover rate in the accumulation period first decreased, then increased. In the stable period, the snow cover rate was stable at 15% to 20%. During the ablation period, the snow cover rate first decreased, then increased and then decreased. The spatial distribution of snow cover in the Babao River Basin varies greatly, and the spatial distribution of snow cover in snow cover periods is different. During the accumulation period, the snow cover is only distributed in the mountainous areas. During the stable period, the snow cover area does not change much. The snow cover is stable in the high-altitude mountains and river valleys, the snow changes rapidly.