Semi-Automated Mapping of Complex-Terrain Mountain Glaciers by Integrating L-Band SAR Amplitude and Interferometric Coherence

: Mapping the outlines of glaciers has primarily relied on the interpretation of satellite optical images. However, the accurate delineation of glaciers in complex terrain mountain regions remains challenging, mainly because the supraglacial debris-covered ablation zones and snow-covered accumulation zones often exhibit the same spectral properties as their adjacent grounds in optical images. This study presents a novel approach by exploring both the satellite synthetic aperture radar (SAR) amplitude and interferometric coherence to map mountain glaciers. This method explores the deviation of the glacier surface signal in the SAR time series to distinguish glacier ice from the surrounding stable ground. To this end, we explored the classifying capabilities of two indices from a set of SAR images, SAR interferometric coherence and amplitude deviation index (ADI), to determine glacier boundary. We found that the two indices complement each other for mapping glaciers. A ratio map based on ADI and SAR coherence (ACR) was then derived, from which the glacier outline was automatically tracked using a speciﬁed threshold, followed by manual modiﬁcation. We validated this approach on two typical valley glaciers, the debris-covered Hailuogou Glacier and debris-free Mozigou Glacier, in Mount Gongga in the southeastern Tibetan Plateau. The results show that the proposed ACR criteria can signiﬁcantly enhance the contrast between glaciers and their surroundings. By comparing our results with manually delineated glacier outlines from high-resolution cloud-free satellite optical imagery, we found that the misclassiﬁcation rate and difference rate for our results were 2.6% and 4.2%, respectively. The approach presented in this study can be easily adapted to map the outlines of mountain glaciers worldwide efﬁciently and is useful for inferring glacier boundary changes in a climate warming context. on various parts of the glacier. selected two pure-ice glacier patches in the accumulation region (GA1 and two bedrock patches around the accumulation region (BA1 and BA2), two debris-covered glacier tongue region (GT1 GT2), two soil near the tongue region


Introduction
Glaciers are an important source of freshwater on earth and an indicator of climate change. Several global glacier inventory initiatives have been launched in recent decades, such as the Global Land Ice Measurements from Space (GLIMS), the Randolph Glacier Inventory (RGI), and the Glacier Area Mapping for Discharge in Asian Mountains (GAM-DAM) Glacier Inventory (GGI). These global glacier inventories provide accurate and Figure 1. Topographic map (a) and Google Earth optical image (b) of the test area. The white (a) and yellow (b) lines illustrate the glacier boundaries from SCGI and GGI, respectively. (c,d) photographs of HLG-1 taken at the stations shown as "S1" and "S2" of (b), respectively.

Data Sets
We used eight L-band ALOS-2 PALSAR-2 images, captured between 1 December 2017 and 21 September 2018, to estimate the coherence and ADI. All images were acquired in Stripmap mode 1 (SM1) and HH polarization from an ascending path (Path 147), with a time span of 42 days for each pair of adjacent images. The pixel spacing for the SAR images is 1.83 × 1.43 m in the SAR azimuth and range directions, respectively. We selected the Lband SAR data because of its stronger penetrability, which has the advantage of overcoming phase decorrelation caused by seasonal changes in vegetation-covered areas [29].
SAR data processing requires an external digital elevation model (DEM) to simulate and remove the topographic phase. We used the Shuttle Radar Topography Mission (SRTM) DEM with a spatial resolution of approximately 30 m to assist InSAR processing. The DEM was also used as a reference to achieve the transformation from range-Doppler coordinates to the geospatial system. We also collected a high-resolution orthophoto (~10 cm resolution) covering the HLG glacier tongue and a partial icefall area as a base map to verify the accuracy of the glacier boundaries. The orthophoto was produced from DJI drone images collected on 5 June 2019.
We downloaded all cloud-free Sentinel-2 and Landsat-8 satellite optical remote sensing images acquired in 2017 and 2018. However, there is no available data during the warm season. For these cold season images, the Sentinel-2 satellite image acquired on 6 November 2017, which we selected for comparison and discussion, is the least affected by snow cover. In addition, high-resolution Google Earth images captured on 3 December 2011, 29 November 2014, and 28 November 2017, were used to assist in the delineation of the glacier outline. For validation purposes, we used outlines from the SCGI [4] and the GGI [5,6]. Figure 2 shows the data processing flowchart for the proposed method. The main steps of the proposed approach are the co-registration of SAR images, generation of the average coherence map and ADI map, calculation of the ADI-coherence ratio (ACR) map, and delineation of the glaciers. In the following text, we will first describe the principles of using SAR coherence (Section 3.1) and ADI (Section 3.2) to classify glacier. Considering the advantages and disadvantages of these two methods, we will then describe the ACR map as a new criterion for glacier classification. Finally, we show the use of an improved threshold segmentation method to automatically classify glaciers based on the ACR map.

Coherence
SAR coherence represents a statistical criterion for the interferometric phase quality of an SAR image pair [30]. The coherence coefficient is usually estimated by spatially averaging the radar echoes in a moving window [31]: where γ denotes the coherence coefficient of the pixel at coordinates (m, n); µ 1 and µ 2 denote the primary and co-registered secondary images, respectively; M and N define a multi-looking window; and * denotes the complex conjugate operator.

Coherence
SAR coherence represents a statistical criterion for the interferometri of an SAR image pair [30]. The coherence coefficient is usually estimated eraging the radar echoes in a moving window [31]: where γ denotes the coherence coefficient of the pixel at coordinates (m, n note the primary and co-registered secondary images, respectively; M multi-looking window; and * denotes the complex conjugate operator. The value of SAR coherence is influenced by several factors [31], su decorrelation, temporal decorrelation, volume decorrelation, thermal no noise, and atmospheric delay. Assuming that the thermal noise, atmosph ionospheric noise are ignored, a coherence value of 1 indicates no change tion of scatterers within the estimation window during the image acquisit 0 represents complete decorrelation. Owing to the rapid movement and glacier, the coherence coefficient corresponding to the glacier is generally The value of SAR coherence is influenced by several factors [31], such as geometric decorrelation, temporal decorrelation, volume decorrelation, thermal noise, ionospheric noise, and atmospheric delay. Assuming that the thermal noise, atmospheric delay, and ionospheric noise are ignored, a coherence value of 1 indicates no change in the composition of scatterers within the estimation window during the image acquisitions. The value 0 represents complete decorrelation. Owing to the rapid movement and melting of the glacier, the coherence coefficient corresponding to the glacier is generally severely decorrelated. In contrast, the coherence coefficient is generally higher for relatively stable non-glacial objects. We used the commercial GAMMA software to process the SAR images [32]. Based on the eight SAR images, we obtained 12 image pairs with a maximum time interval of 84 days and a perpendicular baseline limit of 240 m. The details of the selected InSAR pairs are shown in Table 1. We applied the multi-looking operation with six looks in the range direction and eight looks in the azimuth direction and filtered each interferogram using the adaptive Goldstein filter with a window size of 8 × 8 pixels. Finally, the coherence maps of the selected 12 InSAR image pairs calculated using Equation (1) are shown in Figure 3.      As analyzed previously, Figure 3 shows the stable low coherence characteristics of the glacier-covered area, while for non-glacial areas, the coherence has a significant variation in the time series. In addition, we did not observe significant correlations of coherence with temporal and spatial baselines. However, the coherence of image pairs As analyzed previously, Figure 3 shows the stable low coherence characteristics of the glacier-covered area, while for non-glacial areas, the coherence has a significant variation in the time series. In addition, we did not observe significant correlations of coherence with temporal and spatial baselines. However, the coherence of image pairs spanning the cool and warm seasons shows more severe decoherence. This suggests that seasonal ground surface variation may be the most significant factor influencing coherence.
In glacial regions, the coherence would keep low-level coherence due to the largescale displacements. In contrast, the coherence varies significantly in non-glacial regions. Thus, both the average and standard deviation of the coherence maps could be used as an indicator of glacier classification. The average and deviation for the glacier surface are at low levels, while the opposite is true for the non-glacier regions. However, considering that the adjacent permanent scatterers would show a pattern of high mean and low standard deviation, we took an average of the coherences from multiple SAR image pairs to ensure the robustness of the glacier classification criteria.

SAR Amplitude Dispersion Index
Amplitude Dispersion Index (ADI) was presented as a discrimination method for Persistent Scatterer (PS) point selection and applied to Persistent Scatterer Interferometric SAR technology [27], which is defined as: where σ amp and a denote the standard deviation of the amplitude and the mean amplitude, respectively; a i and |a i | represent the complex value and amplitude in the ith image, respectively; and N is the number of images. As a time-domain estimation, the ADI reflects the stability of the backscatter coefficient in the SAR time series. Valley glaciers are generally characterized by instability due to movement and ablation. In contrast to the adjacent environment, the backscattered signal on the glacier surface is extremely unstable, resulting in a higher-level ADI. Under these assumptions, we propose the ADI as an additional criterion to discriminate glaciers from their surroundings.
The amplitude fluctuation of SAR images can be caused by several factors, for example, instability of the ground surface, satellite incident angle, air humidity, thermal noise, and ionospheric noise [33]. In this case, it is assumed that the attenuation of the backscattering intensity of SAR images due to air humidity in a finite space is uniform. In addition, the angle of incidence was fixed for each image. Thus, if the thermal noise and ionospheric noise are neglected, the most significant factor resulting in amplitude fluctuation is the change in the ground surface. Although the humidity and roughness of periglacial objects tend to be more sensitive than glaciers caused by seasonal freeze-thaw and rainfall, the changes in the glacial surface are more dramatic because the material changes play a more prominent role in amplitude fluctuation.

Generation of an ACR Map Combining ADI and Coherence
Using Equation (1), we calculated the coherence coefficient for each pixel based on a window size of 7 × 7 pixels. To calculate the ADI, we first estimated the temporally averaged amplitude and standard deviation of each pixel, and then obtained the ADI map using Equation (2). Figure 4a,b show the averaged coherence and ADI maps, respectively. Both visually exhibit the capability of glacier classification.
Taking the HLG-1 as an example, the distributions of low coherence and moderate/high ADI maps are remarkably correlated for the majority of glacier extent. However, we observe that the coherence of glaciers and non-glaciers in the accumulation zone is at the same level. Therefore, it is difficult to directly achieve the identification of glaciers solely using the coherence map. In comparison, the contrast between the glacial ice and exposed bedrock is extremely obvious on the ADI map. Conversely, for the glacier ablation region (marked on Figure 3 as the glacial tongue) the ADI values of glacial and non-glacial areas are almost equal, while the coherence map shows the advantages of glacier classification.
To quantitatively express the ability of coherence and ADI maps to classify glaciers in accumulation and ablation regions, we selected eight equally sized rectangular patches (51 × 21 pixels) based on the GGI outline ( Figure 4a) in and around HLG-1. The sample patches were located on various parts of the glacier. We selected two pure-ice glacier patches in the accumulation region (GA1 and GA2), two bedrock patches around the accumulation region (BA1 and BA2), two debris-covered glacier patches in the tongue region (GT1 and GT2), and two bare soil patches near the tongue region (BT1 and BT2). Figure 4c,d shows the two-dimensional distribution maps of the ADI and coherence over the eight patches. The sample patches corresponding to glaciers and non-glaciers were clustered. In the HLG-1 accumulation region, the coherence range of glacier and nonglacier partially overlap, while their ADI distributions are entirely separated. However, this pattern was the opposite in the ablation region. There is an apparent intersection in the ADI value range of the samples, while the glacier and its surroundings are clearly divided in the coherence range.
tion region (marked on Figure 3 as the glacial tongue) the ADI values of glacial and nonglacial areas are almost equal, while the coherence map shows the advantages of glacier classification.
To quantitatively express the ability of coherence and ADI maps to classify glaciers in accumulation and ablation regions, we selected eight equally sized rectangular patches (51 × 21 pixels) based on the GGI outline ( Figure 4a) in and around HLG-1. The sample patches were located on various parts of the glacier. We selected two pure-ice glacier patches in the accumulation region (GA1 and GA2), two bedrock patches around the accumulation region (BA1 and BA2), two debris-covered glacier patches in the tongue region (GT1 and GT2), and two bare soil patches near the tongue region (BT1 and BT2). Figure 4c,d shows the two-dimensional distribution maps of the ADI and coherence over the eight patches. The sample patches corresponding to glaciers and non-glaciers were clustered. In the HLG-1 accumulation region, the coherence range of glacier and non-glacier partially overlap, while their ADI distributions are entirely separated. However, this pattern was the opposite in the ablation region. There is an apparent intersection in the ADI value range of the samples, while the glacier and its surroundings are clearly divided in the coherence range.  From the above analysis, we found significant patterns of high-level ADI and low-level coherence for the glaciers, while the opposite was true for the non-glacial objects. Thus, utilizing the complementary and reverse features of coherence and ADI, we performed pixel-by-pixel ratio operations with ADI as the numerator and coherence as the denominator to generate the ADI-coherence ratio (ACR) index (Equation (3)). The calculated ACR map is shown in Figure 5a. The resolution of the ACR map is consistent with the coherence and the ADI maps and is related to the raw resolution and multi-look factor of the SAR image (about 11 m in this case). Compared with the ADI and coherence maps, the ACR map maximizes the contrast between the glacier and its surroundings. It integrates the advantages of coherence and ADI to distinguish the glacier from the glaciation region, regardless of the accumulation and ablation portions of the glacier.
where ADI(r, c) and Coh(r, c) denote the ADI value and coherence of the pixel (r, c) in the image space coordinate system, respectively.
denominator to generate the ADI-coherence ratio (ACR) index (Equation (3)). The calculated ACR map is shown in Figure 5a. The resolution of the ACR map is consistent with the coherence and the ADI maps and is related to the raw resolution and multi-look factor of the SAR image (about 11 m in this case). Compared with the ADI and coherence maps, the ACR map maximizes the contrast between the glacier and its surroundings. It integrates the advantages of coherence and ADI to distinguish the glacier from the glaciation region, regardless of the accumulation and ablation portions of the glacier.

ADI( , ) ACR( , ) Coh( , )
r c r c r c = (3) where ADI(r,c) and Coh(r,c) denote the ADI value and coherence of the pixel (r,c) in the image space coordinate system, respectively.
To verify the effectiveness of the ACR map for distinguishing the glacier from its surroundings, we used histograms to count the ACR values of the exemplified patches marked in Figure 4a. Figure 5c shows the histograms for the four patches in and around the accumulation region, and Figure 5d shows the histograms for the glacier tongue region. The statistical histograms of these patches illustrate the applicability of the ACR as a glacier classification criterion. Especially in the accumulation region, the sample patches in the glacier and non-glacial areas were completely partitioned to different ACR levels. For the glacial tongue, patches from glacial and non-glacial samples were also effectively clustered, despite the ACR value range of these sample patches being fractionally overlapped. Table 2 shows the statistics and comparisons of ADI, coherence and ACR attributes for the selected samples.  To verify the effectiveness of the ACR map for distinguishing the glacier from its surroundings, we used histograms to count the ACR values of the exemplified patches marked in Figure 4a. Figure 5c shows the histograms for the four patches in and around the accumulation region, and Figure 5d shows the histograms for the glacier tongue region. The statistical histograms of these patches illustrate the applicability of the ACR as a glacier classification criterion. Especially in the accumulation region, the sample patches in the glacier and non-glacial areas were completely partitioned to different ACR levels. For the glacial tongue, patches from glacial and non-glacial samples were also effectively clustered, despite the ACR value range of these sample patches being fractionally overlapped. Table 2 shows the statistics and comparisons of ADI, coherence and ACR attributes for the selected samples. The selection of a threshold is a critical step for delineating the outlines of glaciers based on the ACR map. In the above section, we showed that the ACR map derived from multitemporal SAR images can be used to classify glaciers. To achieve semi-automatic delineation of glacier boundaries, we proposed a process of automatic binarization classification based on the ACR map and then manually mapped glacier outlines.
The ACR map is generated using the ratio of the two arrays, and the distribution of the ACR values is close to the exponential function. We applied the logarithmic transformation of the value of the ACR map to compress its value range and then performed a normalization operation. Figure 5b shows a histogram of the log-transformed ACR map, which exhibits an apparent bimodal distribution. Data with a bimodal distribution can be identified as an ensemble of two separable groups. However, glaciers correspond to different ACR values. Even local points on the same glacier vary in ACR values, depending on various parameters such as the slope and mass. Figure 5c,d show that the range of ACR values is different between the glacial accumulation and ablation regions. A fixed global threshold for classifying the ACR map would lead to a misclassification of glaciers. Thus, we introduced an adaptive threshold segmentation method to improve this process [34].
To reduce the misclassification of glaciers, we propose a method that combines global thresholding and local thresholding to automatically achieve segmentation of the ACR map. We first performed Otsu's method to generate the threshold T otsu for separating the ACR map into two glacier and non-glacier classes [35]. Otsu's method can automatically find the optimal threshold to maximize the variance between classes and has been widely applied for image binarization processing. We then filtered the normalized ACR map with 199 × 199 windows using an average function and multiplied the filtered ACR map by 0.9 as the adaptive threshold matrix T local . For the coordinate (r, c) of the ACR map, the final threshold T(r, c) is defined as: where T local (r, c) and T otsu denote the thresholds derived from the adaptive method and Otsu's method, respectively; r and c denote the row and column numbers of pixels in the SAR image space coordinate system, respectively. After the above steps, we converted ACR to a binary map based on the generated threshold array and removed all connected objects with fewer than 99 pixels from the binarized ACR map. We used the DEM to achieve the projection transformation from the range-Doppler coordinate to the geographic coordinate for the binarized ACR map, as shown in Figure 6a. Considering the inherent geometric distortion of the SAR image, we calculated the layover-shadow map using the DEM and satellite parameters, as shown in Figure 6b.
Accurate automatic mapping of glacier boundaries is difficult because of the rough edges of the binarized ACR maps. Thus, we manually delineated and refined the glacier boundaries by combining a binarized ACR map and high-resolution satellite optical images. We first filtered the non-glacial raster pixels in the binarized ACR map, masked the areas corresponding to layover and shadow, and then overlaid these layers on the high-resolution Google Earth satellite image (Figure 6c). Although glaciers and surroundings were classified in the binarized ACR map, multiple dendritic branches adjacent to glaciers caused by unstable slopes needed to be filtered. Thus, we manually traced the outline of the glacier along the edge of the classified ACR map and delineated the exposed bedrock in the interior of the glacier. Identification of the features of each classified object, such as hanging glaciers and avalanche chutes, was verified visually using Google Earth images. The boundaries of adjacent glaciers were determined manually by topography on Google Earth.
The initial mapping of the glacier boundary must be followed by a refinement operation. Data-free areas affected by layover and shadows were identified using optical images. Since most of these disturbed areas exist in the glacial accumulation region, we distinguished the classes of these masked surfaces by visual interpretation. Glacier boundaries were corrected accordingly. Finally, the attributes of each glacier, such as terminal elevation and area, were calculated with the aid of the GIS tool. resolution Google Earth satellite image (Figure 6c). Although glaciers and surroundings were classified in the binarized ACR map, multiple dendritic branches adjacent to glaciers caused by unstable slopes needed to be filtered. Thus, we manually traced the outline of the glacier along the edge of the classified ACR map and delineated the exposed bedrock in the interior of the glacier. Identification of the features of each classified object, such as hanging glaciers and avalanche chutes, was verified visually using Google Earth images. The boundaries of adjacent glaciers were determined manually by topography on Google Earth. The initial mapping of the glacier boundary must be followed by a refinement operation. Data-free areas affected by layover and shadows were identified using optical images. Since most of these disturbed areas exist in the glacial accumulation region, we distinguished the classes of these masked surfaces by visual interpretation. Glacier boundaries were corrected accordingly. Finally, the attributes of each glacier, such as terminal elevation and area, were calculated with the aid of the GIS tool.

Accuracy Assessment
To validate the effectiveness of the proposed method, we manually mapped the glacier boundaries as ground truth using high-resolution Google Earth historical images. With the benefits of high-resolution 3D visualization of Google Earth, we can accurately track the glacier outlines by interpreting the texture differences between the glacier and adjacent bedrock. We mapped glacier boundaries primarily based on Google Earth images acquired in November 2017 to attenuate the potential bias caused by glacier changes. Considering that the median date of SAR image acquisitions is 27 April 2018, we assume that the changes in glacier area during this 150-day interval are ignored. We then used Google Earth historical images to cross-identify the exposed bedrock and modify the extent of the glacier. We followed the GLIMS Analysis Tutorial and filtered out snow-adhered rock walls above the bergschrund in the accumulation region. All mapping of glacier boundaries was performed in the WGS-84 geographic coordinate system reduce the effect of projection errors between the heterogenous remote sensing imagery involved. The area counts and comparisons were uniformly projected onto the UTM 47N coordinate system for calculation.
We present the area difference rate, misclassification rate, and deficiency rate as the criteria to define the accuracy of the boundaries derived from the proposed method and the ground truth. The difference rate is the relative area error between the item to be evaluated and the ground truth, calculated by dividing the absolute error by the magnitude

Accuracy Assessment
To validate the effectiveness of the proposed method, we manually mapped the glacier boundaries as ground truth using high-resolution Google Earth historical images. With the benefits of high-resolution 3D visualization of Google Earth, we can accurately track the glacier outlines by interpreting the texture differences between the glacier and adjacent bedrock. We mapped glacier boundaries primarily based on Google Earth images acquired in November 2017 to attenuate the potential bias caused by glacier changes. Considering that the median date of SAR image acquisitions is 27 April 2018, we assume that the changes in glacier area during this 150-day interval are ignored. We then used Google Earth historical images to cross-identify the exposed bedrock and modify the extent of the glacier. We followed the GLIMS Analysis Tutorial and filtered out snow-adhered rock walls above the bergschrund in the accumulation region. All mapping of glacier boundaries was performed in the WGS-84 geographic coordinate system reduce the effect of projection errors between the heterogenous remote sensing imagery involved. The area counts and comparisons were uniformly projected onto the UTM 47N coordinate system for calculation.
We present the area difference rate, misclassification rate, and deficiency rate as the criteria to define the accuracy of the boundaries derived from the proposed method and the ground truth. The difference rate is the relative area error between the item to be evaluated and the ground truth, calculated by dividing the absolute error by the magnitude of the ground truth area. The misclassification rate and deficiency rate are defined as calculating the spatial relation between the ACR-based glacier boundaries and ground truth. For outlines A and B, we divide the area into four categories, not in A nor in B (C1), only in A(C2), only in B (C3), and both in A and B(C4), as shown in Equation (5), and calculate the percentage of their area.
where C1 to C4 indicate sub-categories; p denotes any element of the set containing the subset A and B; "&" denotes "and".
We also compared the boundaries from the SCGI and GGI with the ground truth and calculated their areal difference rate, misclassification rate, and deficiency rate. Then, the accuracy corresponding to SCGI and GGI was compared horizontally with that of the proposed method. Figure 7 shows the interpreted boundaries following the method described above, in which the base map is a Sentinel-2 false-color image acquired on 6 November 2017. We identified all glaciers in the HLG and MZG catchments (i.e., nine and four, respectively). Some internal polygons represent exposed bedrock that exclude glacial extent. Glacier details are shown in Table 3.

Validations of the Glacier Outlines from the ACR Map
We manually mapped the boundary of HLG-2 as the ground truth based on Google Earth historical images acquired in December 2011, November 2014, and November 2017. We selected HLG-2 as the validation area because it is less affected by the SAR imaging geometric distortions such as layover and shadow (Figure 6b). HLG-2 is a debris-covered glacier with a large amount of ice mass detached from the headwall in the accumulation region and several exposed bedrock areas in its interior. We visually interpreted the nonglacial areas and manually traced the glacial outlines based on high-resolution images (Figure 8). The bare ground in the interior of the glacier was also identified and delineated to ensure delineation accuracy (yellow polygons in Figure 8a). Figure 8b,c show the true outlines of HLG-2 superimposed on the ACR map and the binarized ACR map, respectively. Figure 8d shows the glacier boundaries from four datasets (GGI, SCGI, ACR, and the true data) overlaid on the Google Earth image. In visual

Validations of the Glacier Outlines from the ACR Map
We manually mapped the boundary of HLG-2 as the ground truth based on Google Earth historical images acquired in December 2011, November 2014, and November 2017. We selected HLG-2 as the validation area because it is less affected by the SAR imaging geometric distortions such as layover and shadow (Figure 6b). HLG-2 is a debris-covered glacier with a large amount of ice mass detached from the headwall in the accumulation region and several exposed bedrock areas in its interior. We visually interpreted the nonglacial areas and manually traced the glacial outlines based on high-resolution images (Figure 8). The bare ground in the interior of the glacier was also identified and delineated to ensure delineation accuracy (yellow polygons in Figure 8a). Figure 8b,c show the true outlines of HLG-2 superimposed on the ACR map and the binarized ACR map, respectively. Figure 8d shows the glacier boundaries from four datasets (GGI, SCGI, ACR, and the true data) overlaid on the Google Earth image. In visual comparison with the ground truth, the outline from our ACR map exhibits the best agreement with the true data compared with the GGI and SCGI datasets. SCGI had the most notable localization errors. To quantitatively evaluate the accuracy of the proposed method, we further calculated the deviation, difference rate, misclassification rate, and deficiency rate of the boundaries derived from the SCGI, GGI, and ACR maps with respect to the ground truth ( Table 4). The spatial relationships between ground truth and glacier boundaries of our result, GGI, SCGI are shown in Figure 9.
The area determined from the ACR map was smaller than 11% and 14% of those from the GGI and SCGI, respectively. This indicates that the GGI and SCGI may have overestimated the glacier-covered area for Mount Gongga. The area difference of glaciers generated from the ACR map and ground truth was the smallest, with a difference rate of 4.4%. The misclassification rate and difference rate for the ACR boundary are 2.6% and 4.2%, respectively, which was close to those for the GGI of 3.9% and 2.7%, respectively. This is likely owing to localization bias induced by the SAR image projection transformation. In addition, note that the accuracy of glacier boundaries from the SCGI in this region is significantly lower than that of the GGI and ACR results. comparison with the ground truth, the outline from our ACR map exhibits the best agreement with the true data compared with the GGI and SCGI datasets. SCGI had the most notable localization errors. To quantitatively evaluate the accuracy of the proposed method, we further calculated the deviation, difference rate, misclassification rate, and deficiency rate of the boundaries derived from the SCGI, GGI, and ACR maps with respect to the ground truth ( Table 4). The spatial relationships between ground truth and glacier boundaries of our result, GGI, SCGI are shown in Figure 9. The area determined from the ACR map was smaller than 11% and 14% of those from the GGI and SCGI, respectively. This indicates that the GGI and SCGI may have overestimated the glacier-covered area for Mount Gongga. The area difference of glaciers generated from the ACR map and ground truth was the smallest, with a difference rate of 4.4%. The misclassification rate and difference rate for the ACR boundary are 2.6% and 4.2%, respectively, which was close to those for the GGI of 3.9% and 2.7%, respectively. This is likely owing to localization bias induced by the SAR image projection transformation. In addition, note that the accuracy of glacier boundaries from the SCGI in this region is significantly lower than that of the GGI and ACR results.

Discussion
Our test on the HLG and MZG showed the capability of the ACR model to enhance the contrast between the glacier and its surroundings. Next, we will discuss the comparison of ACR with optical remote sensing for interpreting glaciers and the advantages of ACR over solely using coherence. Finally, we analyze and summarize the limitations of the ACR model for glacier delineation.

Comparisons with the Glacier Classifications from Satellite Optical Images
Although vulnerable to cloud and daylight shadows, satellite optical images have been the primary data sources for glacier mapping over the last decades [3][4][5][6]. In the SWIR band, the reflectance of snow and ice is very low, which is widely exploited for snow classification. The reflectivity of a false-color image using the SWIR-2, green, and red bands is a good indicator of ice and bare ground, which takes advantage of the SWIR enhancement of snow and ice [36,37], as shown in Figure 7. Pure ice and snow show as a cyan color, whereas bare ground shows as a red color. In the glacier accumulation region, the ACR-derived glacier boundaries are consistent with the ice and snow cover boundaries in the false-color image. However, the irregular shape of the outlines in the accumulation region reflects the complex ice distribution. The presence of seasonal snow in the accumulation region probably blurs the ice-rock boundary in the false-color image because the Sentinel-2 optical satellite image was acquired in the cold season. In the glacial tongue zone of HLG-1, the glacier boundary cannot be distinguished based merely on the SWIR band data due to the surface debris cover.
We compared the ACR-derived glacier outlines with the present glacier inventories GGI and SCGI, both derived from the visual interpretation of optical satellite images [4,5]. Figure 10a shows the glacier boundaries from the SCGI, GGI, and this study overlaid on the Google Earth image. The glacier boundaries produced by the different methods were generally consistent. However, noticeable disparities were also visible in some parts of the

Discussion
Our test on the HLG and MZG showed the capability of the ACR model to enhance the contrast between the glacier and its surroundings. Next, we will discuss the comparison of ACR with optical remote sensing for interpreting glaciers and the advantages of ACR over solely using coherence. Finally, we analyze and summarize the limitations of the ACR model for glacier delineation.

Comparisons with the Glacier Classifications from Satellite Optical Images
Although vulnerable to cloud and daylight shadows, satellite optical images have been the primary data sources for glacier mapping over the last decades [3][4][5][6]. In the SWIR band, the reflectance of snow and ice is very low, which is widely exploited for snow classification. The reflectivity of a false-color image using the SWIR-2, green, and red bands is a good indicator of ice and bare ground, which takes advantage of the SWIR enhancement of snow and ice [36,37], as shown in Figure 7. Pure ice and snow show as a cyan color, whereas bare ground shows as a red color. In the glacier accumulation region, the ACR-derived glacier boundaries are consistent with the ice and snow cover boundaries in the false-color image. However, the irregular shape of the outlines in the accumulation region reflects the complex ice distribution. The presence of seasonal snow in the accumulation region probably blurs the ice-rock boundary in the false-color image because the Sentinel-2 optical satellite image was acquired in the cold season. In the glacial tongue zone of HLG-1, the glacier boundary cannot be distinguished based merely on the SWIR band data due to the surface debris cover.
We compared the ACR-derived glacier outlines with the present glacier inventories GGI and SCGI, both derived from the visual interpretation of optical satellite images [4,5]. Figure 10a shows the glacier boundaries from the SCGI, GGI, and this study overlaid on the Google Earth image. The glacier boundaries produced by the different methods were generally consistent. However, noticeable disparities were also visible in some parts of the outline. Compared with GGI and SCGI, our delineations exclude more internally exposed bedrock and bare ground. Consequently, the area of MZG-1 derived from this study was smaller than that of the SCGI and GGI (see Table 2). As an example, Figure 10d shows that more snow-adhered bedrock was defined as non-glacier areas owing to the penetration capability of SAR data. Furthermore, the terminal outline of MZG-1 derived from the ACR map retreated more than 100 m compared with that from the GGI (see Figure 10b), which was delineated based on the Landsat-5 image acquired on 20 August 2000.
A comparable situation was observed for HLG-1. HLG-1 is relatively much more complex because of its glacial tongue covered by moraines and its rugged topography and cirque-type accumulation region. The area of HLG-1 determined from the ACR map is smaller than that derived from SCGI and GGI, because we have identified more nonglacier targets such as internal bedrock and headwalls in the accumulation region, which are excluded from the glacier area calculation. Figure 10e shows glacier outlines in the accumulation region of HLG-1, where a large amount of snow-attached but visibly striated bedrock surface was distinguished as a non-glacial area.
Multispectral-based image enhancement techniques present an irreplaceable advantage in the automated classification of debris-free glaciers. Raup and Khalsa (2010) suggested that snow-covered areas in images acquired at the end of the glacial ablation period could be defined as perennial snow areas and treated as part of the glacier, mainly because multispectral remote sensing images cannot directly interpret ground objects under perennial snow cover [16]. However, optical remote sensing has an inherent flaw of the inability to penetrate snow and ice. In case there is snow attached to the bedrock surface in the accumulation area at high altitude, it will reduce the accuracy of glacier classification and even cause misclassification.
The penetration of SAR images for dry snow makes it possible to distinguish these snow-covered and snow-attached features [38]. As shown in Figure 10e, because most of this region adjacent to the main peak is covered with snow (even the headwalls above the bergschrund have snow-mass adhering to them) it challenging to identify the landforms directly based on spectral feature analysis. Thus, ridgelines are primarily used to determine valley glacier accumulation zone boundaries, such as SCGI. This facilitates the division of independent glaciers but leads to a reduction in the precision of classification, and even to misclassification.
Another advantage of the ACR model over optical remote sensing is that it can avoid daylight shadows. Compared to SCGI, we identified three new hanging glaciers (MZG-2, MZG-3, and MZG-4) on the northern side of MZG-1. Furthermore, the boundaries of HLG-1 and HLG-2 from the SCGI were also observed to be incorrectly mapped. For the glacial tongue of HLG-1, cloud-free satellite images are accompanied by daylight shadow due to its deep valley terrain. This would confuse the manual tracking of moraine-covered glacier boundaries. Both GGI and SCGI have misclassified boundaries in this zone, with the former defining outlines along the ridgeline and the latter having apparent deviations due to shadow influence (Figure 10c).

Comparisons of ACR with SAR Coherence for Glacier Classification
ADI reflects the stability of the backscattering intensity. Glacier movement, melting, and variation in surface water content will cause backscatter intensity changes, resulting in higher dispersion of time-series SAR intensity values. Freezing and melting direct influence on the roughness and humidness of the glacier surface. SAR intensity values on the glacier surface would vary periodically interannually. It is necessary for the time series of SAR images to span both cold and warm seasons when using ADI as an indicator for glacier classification. In addition, considering the inherent speckle noise of SAR, the robustness of ADI is better ensured by as many images as possible. Therefore, we suggest using the maximum number of SAR images in a natural year to calculate the ADI map. Remote Sens. 2022, 14, x FOR PEER REVIEW 17 of 21  Although ADI is much less sensitive to surface changes than SAR coherence, it has the advantage of overcoming spatiotemporal decorrelation. Coherence is more sensitive to changing ground features, whereas ADI is more valid for stable target interpretation. Thus, the ACR model complements the advantages of ADI and coherence and further enhances the contrast between glaciers and non-glaciers.
For the glacial accumulation region, given the penetration of the SAR image and the strong backscattering value of the rock surface, the multi-looked intensity values would maintain a relatively stable level. Thus, the ability of ADI to classify glacier and bare areas for the accumulation region compensates for the shortcomings of the coherence map in this respect. The developed ACR model can accurately identify the exposed internal rocks in the glacier accumulation area and filter out the snow-attached headwalls.
For the ablation region covered by moraines, the dense vegetation growing on both sides of the HLG-1 glacial tongue, paraglacial unstable slopes, and even the smooth rock surface also resulted in severe decorrelation, which certainly increased the difficulty of glacier classification using coherence maps. Moreover, under the influence of a complex and variable climate, the dielectric constant of the ground surface changes accordingly, causing a decrease in coherence or even a complete decorrelation.
The ACR map is sensitive to the changing surface while filtering out objects with a stable backscattering intensity in the glacier-free areas. Taking the tongue of HLG-1 as an example, the radar electromagnetic waves acting on the glacier polish would exhibit specular reflection, which would lead to severe decorrelation and be difficult to distinguish from the glacier in the coherence map. The low but stable intensity features corresponding to these objects would cause a remarkable contrast with the glacier surface in the ACR map.
Most valley glaciers in high mountain Asia are characterized by high altitudes and the presence of vegetative growth at the terminus, such as the eastern sections of the Nyainqentanglha Mountains and Himalayas [39,40]. Dense thickets grow near HLG-1. Thus, SAR coherence-based methods cannot be applied directly. Given the typical representativeness of the glaciers in this study site, the ACR model can provide an objective criterion for inventory and status surveys of cloud-obscured glaciers in the southeastern Tibetan Plateau.

Limitations of Glacier Mapping Based on ACR
The ACR method we propose uses the sensitivities of coherence and ADI to ground surface changes. The tests on Mount Gongga have shown that ACR can be used to classify glaciers as either debris-covered or debris-free glaciers. Nonetheless, the ACR method also has limitations, such as layover and shadow effects, and geometric distortion.
First, geometric distortions induced by the inherent side-looking feature of SAR images result in parts of the area being affected by layover and shadow, especially in areas with severe terrain fluctuations (Figure 6b). Such distortions may cause glacier classification omission. Therefore, the misclassifications induced by these geometric distortions require manual discrimination and modification with the assistance of optical images.
Second, an external DEM is always needed to transform SAR images from the range-Doppler system to WGS-84 geographic coordinates. The registration error and resampling algorithm in this process inevitably introduce systematic deviations and local distortions.
In addition, the ACR model is not a fully automatic method when applied to glacier mapping. The key issue is that unstable slopes within a glacier basin can be seriously misclassified. Indeed, the ACR model is intended to provide a criterion for glacier delineation to avoid subjective mistakes by operators.

Conclusions
In this study, we presented a novel approach to delineate and classify glaciers from their surroundings using satellite SAR images. The SAR coherence and amplitude deviation index of SAR images were integrated as the ACR model to enhance the contrast of glaciers and adjacent objects. We demonstrate that the proposed method is applicable to debris-free glaciers and can be directly applied to that type of surface moraine cover. Given the typical representativeness of the glaciers in this study site, the ACR model can provide an objective criterion for mountain glacier inventory and mapping in the southeastern Tibetan Plateau. By selecting the Hailuogou Glacier and the Mozigou Glacier as two test areas, we mapped the glacier boundaries based on the ACR map from eight L-band SAR images.
We draw the following conclusions. (1) As an alternative criterion that we presented to classify glaciers, ADI can supplement the deficiency of SAR coherence to distinguish glacier ice in the accumulation region. Given the penetration capability of radar waves, ADI has the advantage of sensitive identification of headwalls above the bergschrund and internal bedrock, thus overcoming the weakness of SAR coherence, which is susceptible to decorrelation. (2) The ACR model we developed integrates the advantages of SAR coherence and ADI to distinguish glaciers. The derived ACR map can indicate the coefficient of variation in the glaciation area, providing a criterion for classifying glacial and non-glacial objects. This study shows that the ACR model is applicable in both glacial accumulation and ablation zones, regardless of whether the surface is covered by glacier moraines.
The successful application of the proposed method provides an alternative method for inventorying glaciers for regions where it is difficult to acquire cloud-free optical images. Additionally, we suggest that the proposed method can be further applied to monitor glacier boundary changes in a climate warming context.