Automatic delineation of debris-covered glaciers using InSAR coherence derived from X-, C- and L-band radar data: a case study of Yazgyl Glacier

Abstract

Despite their importance for mass-balance estimates and the progress in techniques based on optical and thermal satellite imagery, the mapping of debris-covered glacier boundaries remains a challenging task. Manual corrections hamper regular updates. In this study, we present an automatic approach to delineate glacier outlines using interferometrically derived synthetic aperture radar (InSAR) coherence, slope and morphological operations. InSAR coherence detects the temporally decorrelated surface (e.g. glacial extent) irrespective of its surface type and separates it from the highly coherent surrounding areas. We tested the impact of different processing settings, for example resolution, coherence window size and topographic phase removal, on the quality of the generated outlines. We found minor influence of the topographic phase, but a combination of strong multi-looking during interferogram generation and additional averaging during coherence estimation strongly deteriorated the coherence at the glacier edges. We analysed the performance of X-, C- and L- band radar data. The C-band Sentinel-1 data outlined the glacier boundary with the least misclassifications and a type II error of 0.47% compared with Global Land Ice Measurements from Space inventory data. Our study shows the potential of the Sentinel-1 mission together with our automatic processing chain to provide regular updates for land-terminating glaciers on a large scale.

1. INTRODUCTION

The accurate knowledge of glacier extent is a prerequisite to many glaciological studies,for example volume change estimates (Kääb and others, 2012; Gardelle and others, 2013) and ice dynamic modeling (Huss and Farinotti, 2012; Fürst and others, 2016). Change rates of glacier front positions are often used variables in studies of climate change (Ranzi and others, 2004; GCOS, 2011) or natural hazards such as glacial lake outburst floods (Bajracharya and Mool, 2009; Hewitt, 2014). Therefore, projects like the Global Land Ice Measurements from Space (GLIMS) or the Randolph Glacier Inventory (RGI), which aim to continuously provide updated glacier outlines in digital vector format, are of high importance (Ranzi and others, 2004; Pfeffer and others, 2014; Paul and others, 2015). However, these global datasets are often not consistent over debris-covered glaciers, as observed in this pilot study over Yazgyl Glacier in the Karakoram (Fig. 1), since they have been generated by different operators or methods. Various techniques exist to map glacier outlines using optical (Paul and others, 2015; Smith and others, 2015) and thermal (Alifu and others, 2015) remote-sensing data. However, thermal data have spatial constraints in order to delineate a debris-covered terminus on a high-resolution scale (Alifu and others, 2015). Optical data are often limited due to the lack of cloud-free acquisitions and good illumination (Smith and others, 2015). More importantly, almost similar spectral signatures of debris-covered ice and surrounding moraines can lead to a time-consuming manual digitalization and debris misclassification (Paul and others, 2004). However, Smith and others (2015) improved the algorithm of Paul and others (2004) and coupled multispectral image classification using Landsat ETM+ and OLI data with elevation, slope and velocity thresholds to map debris-covered glacier outlines. Their reported misclassification ranged between 2% and 10% of the glacier area. Several studies have integrated topographic and thermal data to separate supraglacial and periglacial debris (Bolch and others, 2007; Shukla and others, 2010; Bhambri and others, 2011; Racoviteanu and Williams, 2012). Thermal-based classifications rely on a warmer surface of the periglacial features than the supraglacial surface, which might not be true for a thick supraglacial debris cover due to its insulating effect (Schauwecker and others, 2015). Moreover, classifications which are largely based on topography are limited by the availability of high-resolution DEM.

Fig. 1. Location map of Yazgyl Glacier in high-mountain Asia (Karakoram). The subset on the right shows the lower catchment and its outlines provided by different sources. The white polygon is the buffer zone for error analysis and the red arrows are the flow vectors of the glacier derived from two TerraSAR-X scenes (2009-07-19 and 2009-12-20). (a/b) Source: Esri, NOAA, DigitalGlobe, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN and the GIS User Community. (c) Background image: © CNES/Airbus DS (2008-01-15).

Interferometrically derived synthetic aperture radar (InSAR) coherence (Rodriguez and Martin, 1992; Michel and Rignot, 1999; Li and others, 2001; López-Martinez and Pottíer, 2007) can be used as an alternative to the techniques stated above. InSAR coherence is the complex correlation between two synthetic aperture radar (SAR) images, indicating the (temporal) stability of the backscatter signal. Its numerical value is between 0 and 1, where 0 indicates a completely decorrelated surface and 1 refers to a stable surface. This technique can be used to separate stable and unstable areas, regardless of their surface type. As the technique is based on SAR data, it is almost independent of cloud cover and solar illumination. InSAR coherence from L-band SAR data have been used as a primary estimator to delineate debris-covered glaciers in Alaska (Atwood and others, 2010), the Pamirs (Jiang and others, 2011) and the Himalaya (Strozzi and others, 2010; Frey and others, 2012). Atwood and others (2010) coupled the coherence with slope, size and morphological filters to map glacier outlines. However, coherence images of Himalayan glaciers, used in Frey and others (2012), contained data voids due to steep terrain. Therefore, additional manual editing was required. Wu and others (2012) outlined the challenges in distinguishing the glacier area from the surrounding non-glaciated areas of almost similar coherence in C-band, which they bypassed by means of a further developed texture analysis. Similarly, Robson and others (2015) combined InSAR coherence with object-based image classification (Rastner and others, 2014) to classify debris-covered ice in the western Himalaya, Nepal. Over the years, these above-mentioned techniques have been used to map glacier boundaries, but a precise delineation of debris-covered glacier tongues without manual editing is still not possible for many regions (Frey and others, 2012; Smith and others, 2015). Therefore, the primary objective of this study is to provide an automatic processing chain to delineate debris-covered glaciers. We used InSAR coherence as the primary estimator as well as slope and morphological operations as secondary filters. In a first step, we optimize the processing chain by estimating the influence of parameters such as spatial resolution, coherence window size and topographic phase removal. In a second step, we investigate the potential of X- (TerraSAR-X), C- (Sentinel-1) and L-band (Phased Array L-Band Synthetic Aperture Radar from the Advanced Land Observing Satellite (ALOS PALSAR)) SAR data to map glacier outlines. To assess the accuracy, we compare our results with existing glacier databases.

2. STUDY SITE

We selected Yazgyl Glacier (GLIMS ID: - G075270E36240N) for this study. The glacier is situated in the Karakoram Range, which is a part of the widely spread mountain ranges of Hindu Kush-Karakoram-Himalaya. Yazgyl Glacier (Fig. 1) has a land-terminating stable terminus covering an overall catchment area of 117 km² (Rankl and others, 2014). It originates from an elevation of 7589 m a.s.l. and terminates at 3277 m a.s.l. with an average of 5546 m a.s.l. (Bajracharya and others, 2014). The glacier has a gentle topography with an average slope of 16° while the debris-covered part is almost flat (mean slope 4°) with numerous supraglacial ponds. The outlines of the glacier inventories for Yazgyl Glacier from GLIMS (Bajracharya and others, 2014) and Rankl and others (2014) only partially integrate debris-covered parts into their glacier outlines, whereas the RGI 5.0 glacier outlines (Cogley and others, 2015) also include a significant portion of its northern part (Fig. 1). No change in glacier extent was expected despite different acquisition dates of the radar data (Table 1) due to the stable front position of Yazgyl Glacier as noted by Rankl and others (2014).

Table 1. Specifications of sensors used for estimation of coherence

3. DATA

In order to estimate InSAR coherence, satellite data from three different SAR sensors were chosen, namely the ALOS PALSAR with its Fine Beam (FB) mode, the C-band Sentinel-1 with its Interferometric Wide Swath (IW) mode and the X-band data from TerraSAR-X Strip Map (SM) mode. Details of the datasets are given in Table 1.

The ALOS PALSAR and Sentinel-1 data were acquired with similar temporal baselines while the X-band data used in this study had a relatively large temporal baseline of ~5 months. In general, the longer wavelength of L-band reveals a temporally more stable signal than the short X-band wavelength due to the deeper penetration of the L-band signal into the surface and sensitivity to larger order changes. Longer temporal baselines of InSAR pairs affect the overall image coherence and hence might also impact the co-registration quality. For example, a small amount of solid precipitation or wind drift during the acquisition interval might reduce the co-registration accuracy depending on the radar signal, leading to deteriorated coherence even if there is no surface movement (Weydahl, 2001). In our setup, we did not observe low co-registration accuracy despite the larger temporal baseline in case of X-band data. Nevertheless, it should be considered that better results could be achieved with shorter temporal baselines. In this study, we used radar data from summer time (Table 1) in order to specifically use the temporal decorrelation on the glacier surface (Atwood and others, 2010; Weydahl, 2001). As reference data, we used outlines from the updated global glacier inventories GLIMS (Bajracharya and others, 2014), the RGI version 5.0 (Cogley and others, 2015) and the glacier outlines from Rankl and others (2014) (Table 2). These outlines show differences especially over debris-covered part of Yazgyl glacier which is extensively covered by supraglacial ponds (Fig. 1). Since reference data from GLIMS and Rankl and others (2014) are very similar in the terminus of Yazgyl glacier, we only assessed the accuracy with GLIMS and RGI 5.0 data.

Table 2. Overview of reference data

4. METHODS

4.1. InSAR coherence

InSAR coherence determines the statistical noise of the phase in an interferogram between two SAR images (Michel and Rignot, 1999; Frey and others, 2012). The coherence of two images S ₁ and S ₂ can be calculated by averaging over M × N pixels (coherence window) in the image dimensions m × n, mathematically defined as

(1)$$\vert \gamma \vert = \displaystyle{{\sum\nolimits_{m = 1}^M {\sum\nolimits_{n = 1}^N {S_1(m,n)S_2^* (m,n)} } } \over {\sqrt {\sum\nolimits_{m = 1}^M {\sum\nolimits_{n = 1}^N {{\left \vert {S_1(m,n)} \right \vert }^2} } \sum\nolimits_{m = 1}^M {\sum\nolimits_{n = 1}^N {{\left \vert {S_2(m,n)} \right \vert}^2} } } }}$$

(López-Martinez and Pottíer, 2007). $S_{2}^{\ast}$ is the complex conjugate of S ₂.

If other decorrelation effects like thermal noise, geometric or volume decorrelation (Li and others, 2001; Atwood and others, 2010) are neglected, a coherence value of 1 signifies maximum stability, while 0 shows complete temporal decorrelation of the surface between the two acquisition dates. In the case of the glacier, movement or melt processes cause decorrelation resulting in low coherence, but also supraglacial lakes and ponds can induce local decorrelation. The first two parts of Fig. 2 show the proposed processing chain for estimating coherence from the three different SAR sensors, including the intermediate processing steps (multi-looking factor, co-registration, interferogram formation, optional filtering and topographic phase removal and geocoding). The multi-looking factor (MLF) is the numerical factor by which an input single look complex (SLC) SAR image is averaged in slant-range and azimuth directions. In our processing chain, the resulting SAR amplitude image is used to determine the output resolution of the multi-looked complex interferogram generated from the two co-registered SLC images without further resampling. Before interferogram generation, we co-registered the InSAR image pairs using a cross-correlation optimization technique (Werner and others, 2005). In case of TerraSAR-X and ALOS PALSAR data, with a co-registration accuracy of around 0.2 pixel, the temporal decorrelation should not be biased by more than 5% (GAMMA, 2011). Due to a highly variable doppler centroid throughout the image, a comparable higher co-registration accuracy is required in case of Sentinel-1 IW mode using Terrain Observation with Progressive Scans SAR (TOPSAR). Sentinel-1 IW data consist of three sub-swaths, each of them consisting of several bursts. Therefore, a special procedure for co-registration and a de-bursting step is necessary to join all burst data of Sentinel-1 (Veci, 2015). We processed Sentinel-1 data without multi-looking in open source Toolbox SNAP provided by ESA. After co-registration, the complex interferogram as the pointwise complex multiplication of corresponding pixels in the co-registered InSAR pairs was calculated (Touzi and others, 1999). The subsequent coherence estimation measures the phase noise in the interferogram. Therefore, besides the scattering properties of the target, system and sensor parameters (wavelength, slant range resolution, interferometric baseline, incidence angle, system noise) have an impact on the decorrelation (Frey and others, 2012). To minimize the effect of volume scattering by different surface properties and to extract the temporal change, we used summer scenes where seasonal snow is not present in the ablation area.

Fig. 2. Processing chain to derive glacier outlines using InSAR coherence, slope and morphological operations. The processing chain was implemented in GAMMA (except Sentinel-1 without topographic phase removal in Fig. 4e was performed in SNAP). The post-processing and accuracy assessment for every case was carried out in R. Background image(p7-p9): © CNES/Airbus DS (2008-01-15).

Topographic decorrelation depends on the baseline and the local slope and generates non-stationarity of the phase within the coherence estimation window (Lee and Liu, 2001). Therefore, previous authors suggested subtracting the topographic phase from the complex interferogram in order to account for such effects (Atwood and others, 2010; Jiang and others, 2011). However, we determined the coherence directly from the complex interferogram and in the second step from a differential interferogram where we subtracted the topographic phase simulated by SRTM C-band DEM. Voids in the SRTM DEM were filled with zero in order to have no influence on the topographic phase. Our processing chain was not successful after topographic phase removal at the date of our image processing in SNAP. However, we managed the removal in GAMMA software, but only with a too high MLF for a useful delineation. Nevertheless, the results are shown for every scenario to point out the strong influence of the resolution before coherence estimation. In Table 3, the different combinations of MLF, topographic phase removal and coherence window sizes, presented in this study, are listed. In the first section of results, we present the outlines for TerraSAR-X for different resolutions and coherence window sizes. In the second section, we compare different radar frequencies (X-, C- and L-band) for which we produced coherence images of the same ground resolution (15 m) by varying the MLF to match the output resolution.

Table 3. Values of different parameters for coherence estimation in case of TerraSAR-X (TSX), Sentinel-1 (S-1) and ALOS PALSAR (Palsar) data. The corresponding figure number for every combination is listed

4.2 Methods for post-processing

Following the processing chain proposed by Atwood and others (2010), coherence and slope thresholds were applied. Bolch and others (2007) reported a maximum average slope of 10° for the debris-covered parts of the Mount Everest region in ASTER DEM. Paul and others (2004) applied a maximum gradient of 24° for Oberaletschgletscher in the Swiss Alps. Therefore, we applied a slightly higher threshold of 30° similar to Atwood and others (2010) to hold validity for all regions. Every pixel with a coherence value of <0.2 and a slope value up to 30° was classified as a glacier pixel. In case of Yazgyl Glacier, slope values up to 15° occur even at almost flat (mean slope of ~4°) debris-covered parts, possibly due to ice cliffs. This requires a higher slope threshold than the mean slope value. From the coherence and slope thresholding, we obtained a binary mask consisting of two classes - glacier and no-glacier (p2 in Fig. 2). We subsequently applied binary morphological operations that consist of opening and closing patches based on their sizes (Haralick and Shapiro, 1992). One patch is generally a number of connected pixels of one class (e.g. patch of no-glacier pixels within the glacier surface and patch of glacier pixels outside of glacier surface). For this, we defined the kernel size referring to the number of pixels in both directions in a two-dimensional space. In the closing operation, the patches smaller than the kernel are filled whereas, in the opening operation, patches smaller than the kernel are deleted. In an initial step, we applied a closing operation with a smaller kernel size to include real glacier pixels in close proximity to the main glacier, but to avoid integration of wrong patches (p3 in Fig. 2). In the next step, we performed an opening operation with a larger kernel size to delete wrong patches (p4 in Fig. 2). Finally, the larger kernel size was used for a second closing iteration to fill most of the holes within glacier (p5 in Fig. 2). In one processing path, we masked out layover and shadow areas (p6 in Fig. 2) before converting the raster binary mask into a vector format (p7 in Fig. 2) to give also corrected error values (Fig. 4g). All post-processing steps with associated parameters are illustrated in the third part of Fig. 2.

4.3. Accuracy assessment

We compared our glacier outlines, derived from different sensors and processing parameters, with the reference data. We created a buffer zone of 500 m (white polygon in Fig. 1, p8 in Fig. 2) for the assessment. We distinguish four cases (p9 in Fig. 2): (a) ‘true positive’ refers to the intersection of measured glacier area and the glacier area from the reference dataset (correctly classified). (b) ‘True negative’ refers to pixels that were classified as no-glacier in both the cases (correctly not assigned). (c) ‘False negative’ shows the area classified as glacier only by the reference dataset and (d) ‘false positive’ refers to the glacier area mapped by our processing only. It is quite evident from visual inspection that the reference datasets do not detect the debris-covered part or miss large parts of it. Therefore, the ‘false positive’ area mainly appears in this part of the glacier. We adapted nomenclature from statistical hypothesis testing and defined a type II error which refers to the ratio of ‘false negative’ area and the glacier area from the reference dataset (in contrast to type I error which refers to the ratio of ‘false positive’ area and the glacier area from the reference dataset and is in our case no real ‘error’). When possible, we estimated the type II error twice, firstly considering the whole ‘false negative’ area and secondly the corrected ‘false negative’ after masking out the area influenced by layover and shadow. Areas of layover and shadow show mostly low coherence values, but it is not possible to gain any information about the surface. Therefore, we excluded these areas in the second step of error calculation (named ‘corrected’ error) for a better comparison between the sensors. Our results in Fig. 4g show that in case of Yazgyl Glacier the difference is negligible. As this may not be the case for other study sites and only the error including layover and shadow areas is measuring the overall performance, combining data from different image directions should be considered in this case.

4.4. Velocity measurements

Additionally, we exploited the high-resolution TerraSAR-X scenes for surface velocity measurements. We performed SAR offset tracking over already co-registered InSAR image pair using a cross-correlation optimization technique (Strozzi and others, 2002; Seehaus and others, 2015). A window of 128 × 128 pixels (~251 m × 116 m) patch size with a step size of 25 × 25 pixels (~49 m × 23 m) was used to calculate azimuth and slant range displacements. The erroneous displacement values were discarded using an algorithm defined by Burgess and others (2012). This algorithm discards a displacement vector if its orientation differs by a certain threshold from the average orientation computed over surrounding pixels.

5. RESULTS

5.1. Influence of coherence window size and MLFs (spatial resolutions)

Figure 3 shows the influence of different coherence window sizes (3a–3d) and spatial resolutions (3d–3f) on the glacier outlines derived from TerraSAR-X data. In the first case of adapting the coherence window size (3 × 3 pixels, Fig. 3a), there was very little glacierized area retrieved. With 7 × 7 pixels window size (Fig. 3b), the outlines were improved, but no continuous glacier area could be extracted. In the last two cases (11 × 11 and 15 × 15 pixels, Fig. 3c and 3d), the number of holes were significantly reduced and a more connected glacier outline was formed. In a second step, we kept the coherence window size at 15 × 15 pixels (the best case scenario of coherence window size) and varied the MLF during interferogram generation resulting in a different spatial resolution of the output binary mask. The glacier outlines with 1 × 1 MLF (which means no multi-looking), corresponding to 2 m of output resolution, and 15 × 15 MLF, corresponding to 30 m output resolution, are shown in Fig. 3e and 3f, respectively. Fig. 3g and 3h show the type II error comparing our results with the RGI 5.0 and GLIMS outlines for the different combinations of coherence window size and multi-looking. Type II errors are the smallest in case of no multi-looking (2 m resolution). Larger coherence window sizes improved the accuracy (Fig. 3g and 3h) even for high MLF corresponding to 15 m and 30 m of spatial resolution.

Fig. 3. (a)–(f) Glacier outlines derived with different coherence window sizes (3, 7, 11 and 15) and MLFs (1 × 1, 8 × 7, 15 × 15) from TerraSAR-X data. (g) Plot showing the Type II error (%) (corrected by areas of layover and shadow) for different coherence window sizes and MLFs (spatial resolution) with RGI 5.0 as a reference and (h) GLIMS as a reference. Background image: © CNES/Airbus DS (2008-01-15).

Fig. 4. (a/d): Glacier outlines derived from TerraSAR-X, (b/e): Sentinel-1 and (c/f): ALOS PALSAR data for a 15 m output resolution, except 20 m of Sentinel-1 in (b). In (a)–(c), the topographic phase was subtracted from the complex interferogram before coherence estimation, whereas in (d)–(f) no topographic phase removal (without TPR) was applied. (g): Plot showing the type II error (%) (with and without correction by layover and shadow areas) in comparison with RGI 5.0 and GLIMS as a reference. (h): Plot showing the amount of area and corresponding proportion categorized in four different classes (false positive, false negative, true positive and true negative) in comparison with RGI 5.0 and GLIMS as a reference outlines. Background image: © CNES/Airbus DS (2008-01-15).

5.2. Influence of radar frequency and topographic phase

In order to assess the potential of X-, C- and L-band radar data to derive glacier outlines, we adapted the MLFs for different sensors to match an output resolution of 15 m. The coherence window size was kept the same (15 × 15 pixels). We also determined the influence of topographic phase removal (Fig. 4). Keeping the topographic phase in our processing chain does not strongly influence the type II error for Yazgyl Glacier as can be seen for TerraSAR-X and ALOS PALSAR (Fig. 4g). However, TerraSAR-X derived results after topographic phase removal showed some elongated disturbances in the south-eastern part of the glacier. The outlines from TerraSAR-X cover most of the northern debris-covered parts, but outlines from ALOS PALSAR terminate earlier. In contrast, outlines from ALOS PALSAR show a wider terminus, whereas outlines from TerraSAR-X miss even some clean-ice parts of the glacier. The results from Sentinel-1 processed in SNAP software revealed promising results even without removing the topographic phase (Fig. 4e). A very low type II error (3.92% with RGI 5.0 as reference and 0.47% with GLIMS as reference, corrected by layover and shadow, Fig. 4g) was obtained. As we applied processing of Sentinel-1 data in GAMMA only with a too high MLF, this led to an overall increase of type II error of more than 20% while comparing with both reference datasets (Fig. 4g). Considering RGI 5.0 as reference, only both TerraSAR-X products and the Sentinel-1 outlines without topographic phase removal show a higher amount of ‘false positives’, which is an indicator for outlining the debris-covered part farther than RGI 5.0. The overall performance of TerraSAR-X was improved at a higher spatial resolution of 2 m (Fig. 3e), resulting in a type II error (corrected by layover and shadow) of 3.63% with RGI 5.0 and 0.88% with GLIMS as reference (Fig. 3g and 3h). In this setup, the amount of ‘false positives’ over the debris-covered part is comparable with the results from Sentinel-1 in Fig. 4e, even if there are some minor misclassified areas in case of TerraSAR-X.

6. DISCUSSION

6.1. Influence of coherence window sizes and MLFs (spatial resolutions)

Our results show that the spatial resolution and the coherence window size strongly influence the quality of the glacier outlines. Both the generation of the multi-looked interferogram and the subsequent coherence estimation have an impact on the spatial accuracy. Even if the neighboring pixels are not combined into one pixel, the coherence estimation calculates the average of all neighboring pixels in a coherence window which leads to higher uncertainties especially in the border areas of the glacier where the coherence values are a mixture of stable terrain and glacier. Resampling to a higher resolution during post-processing would not regain the prior information. The averaging only led to a high error in coherence estimation if it is based on a lower resolution image (8 × 7 MLF and 15 × 15 MLF, Fig. 3g). On the contrary, in case of no multi-looking (1 × 1 MLF), type II error remains stable on a low value also for larger coherence window sizes. This is in good agreement with results from Wegmüller and Werner (1995) who recommend larger window sizes for areas with medium or low interferometric correlations. In this case, a high standard deviation in a small averaging window would cause a bias to a higher coherence. In our case, the standard deviation decreased from 0.21 for 3 × 3 pixels coherence window size to 0.17 for 15 × 15 pixels within the glacier outlines. Simultaneously, the mean value decreased from 0.45 to 0.18. Also, other error sources affecting non-stationarity of the phase can be minimized by larger window sizes. Nevertheless, as our processing chain has only to distinguish between the two classes on- and off-glacier, and for example the impact of residual topographic phase on coherence bias increases with larger window sizes, we do not recommend larger window sizes to ensure accuracy at the glacier borders.

6.2. Influence of radar frequency and topographic phase

In our inter-comparison of glacier outlines derived from X-, C- and L-band SAR data no direct dependency on radar frequency alone was recognizable. Due to the longer wavelength, L-band data have a much higher penetration depth and therefore a higher proportion of volume scattering. Nevertheless, Huang and others (2017) measured that the portion of surface scattering (56%) on debris-covered ice for ALOS PALSAR is still slightly higher than volume scattering (40%). In the case of the C-band, they found surface scattering of ~ 70% of the total power. Therefore, it is difficult to find correlation as other parameters like the occurrence of water-filled ponds and the following decorrelation or the debris thickness have a strong impact. Moreover, the temporal baselines, but also the initial spatial resolutions, are different between the sensors.

Results from Sentinel-1 processed in GAMMA (Fig. 4b) are not fully suited for inter-comparison due to processing with only a very high MLF, which caused the glacier frontal part to be detected as eroded and disjointed. The type II error in case of C-band Sentinel-1 data with 15 m resolution processed in SNAP (Fig. 4e) is comparable with the type II error of the best possible resolution of TerraSAR-X data with 2 m (Fig. 3e). The amount of misclassified holes in the final glacier outlines is lower in case of Sentinel-1, which indicates a comparably stable and smooth coherence estimation over the glacier surfaces for C-band radar data. C-band has an intermediate of X- and L-band radar frequencies and probably the best combination of radar signal sensitivity to the surface (highest for X-band) and penetration depth (highest for L-band). The ‘false positive’ area derived from the glacier outlines in Fig. 4 shows smaller extents at the debris-covered part for ALOS PALSAR and TerraSAR-X (1.04 km² in Fig. 4f and 0.88 km² in Fig. 4d) than Sentinel-1 (1.69 km², Fig. 4e). However, one has to keep in mind that the area values alone cannot be a measure. Although the ‘false positive’ areas of TerraSAR-X and ALOS PALSAR are close, the glacier outlines are distinct. The glacier outline from TerraSAR-X does not cover the south-eastern terminus but includes the northern part of the terminus. The glacier outline from ALOS PALSAR is comparably wider and includes the south-eastern part but terminates earlier at the northern terminus. In order to compensate for an output resolution of 15 m for inter-comparison, high-resolution TerraSAR-X data had to be strongly multi-looked, leading to dual averaging by multi-looking and coherence estimation. This probably caused the uncertainty at the glacier edges. In contrast, the glacier area over the debris-covered part with TerraSAR-X in 2 m resolution without multi-looking is 2.35 km². The earlier termination of the glacier outline derived from L-band data may be associated with L-band penetration into the glacier surface by mapping lower decorrelation from the subsurface as compared with the surface. Another possibility could be the slow movement at the glacier tongue, which is detected as coherent starting at the longer wavelength of L-band. However, a strict dependency on wavelength is not explained by the results derived from Sentinel-1 in SNAP, as the C-band radar frequency lies in the middle of X- and L-band radar data. The thesis of slow movement is supported by measurements of the glacier surface velocity which reveals a minor flow over the debris-covered part with a magnitude of <0.10 ± 0.0013 m d⁻¹, whereas the upstream glacier flow reaches up to 0.60 m d⁻¹ (Fig. 1). We infer that the debris-covered part of Yazgyl Glacier is almost stagnant and undergoes a down-wasting due to supraglacial ponds, as it is the case for many debris-covered glaciers in high-mountain Asia (Holzer and others, 2015; Vijay and Braun, 2016; Ragettli and others, 2016). This melt also leads to surface decorrelation (low coherence) even if there is no significant horizontal flow. This decorrelation can also appear for other surfaces like snow patches or rock glaciers but is not observed in our study area. Nevertheless, we expect that other interferometric factors (temporal baseline, selection of acquisition date, incidence angle, imaging geometry etc.) are playing a significant role in coherence estimation. Time series analysis in regard to TerraSAR-X coherence has revealed strong seasonal variability. Although the central Karakoram is a relatively dry region, July and August can still experience influence from the summer monsoon.

In our study area, we found a negligible influence of the topographic phase removal for coherence estimation, which may be due to the viewing geometry and gentle topography in our case study. The major discrepancies occur in case of TerraSAR-X data when an elongated south-eastern outline appeared (Fig. 4a) after subtracting the topographic phase before coherence estimation. This is potentially due to an incorrect topographic phase in this area simulated from the SRTM data. Fig. 5 is showing the simulated topographic phase for each of the different sensors. As in each case, one color cycle has a range of 2π, it can be seen that Sentinel-1 with only one cycle (Fig. 5b) within the glacier area would be only slightly affected by a wrong information. In the case of TerraSAR-X (Fig. 5a), the simulation fails in the north-east from the glacier and is probably the explanation of the disturbances at the eastern tongue where an overlay with the glacier is existing. Moreover, the topographic phase simulation for ALOS PALSAR does not appear to correctly reflect the topography. Even if there is no negative impact visible in the final outlines, the coherence images from interferograms with topographic phase subtraction should be treated carefully, especially for larger baselines like in the case of the ALOS PALSAR data. For the case of TerraSAR-X and Sentinel-1 with small perpendicular baselines, we realized that removing the topographic phase is not necessary for glaciers having a flat terminus such as Yazgyl Glacier. In case of different topographic conditions with narrow valleys, the use of high-resolution DEMs for example from TanDEM-X or slope-adaptive range common-band filtering (Santoro and others, 2007) should be considered.

Fig. 5. Simulated topographic phase information displayed in phase cycles of 2π for TerraSAR-X (a), Sentinel-1 (b) and ALOS PALSAR (c) at 15 m output resolution. The topographic phase simulation for ALOS PALSAR does not appear to correctly reflect the topography.

The viewing geometry in the case of ALOS PALSAR and Sentinel-1 is preferable in order to avoid layover and shadow effects. When comparing outlines with similar area, this can be seen in a stronger decrease (6.41% (not corrected) to 4.58% (corrected by layover and shadow)) of type II error for TerraSAR-X in 2 m than for Sentinel-1 in 15 m resolution (4.11% (not corrected) to 3.92% (corrected by layover and shadow), both in comparison with RGI 5.0). The estimated type II errors of best Sentinel-1 and TerraSAR-X products are comparable or even better than results from other studies (0.34–2% in Alifu and others (2015), 2–10% in Smith and others (2015), 11.5% in Rastner and others (2014), 7–17% in Robson and others (2015)). However, this error comparison may not indicate the supremacy of one single algorithm, as most studies, including ours, have been applied to specific test sites only. Moreover, the error assessment depends on the ratio of debris-covered area part to the entire glacier catchment size. In the case of manual digitization by different analysts, a local uncertainty might reach up to 150 m for debris-covered glaciers depending on the available input data (Frey and others, 2012). As a consequence, the proposed automated technique is well within the uncertainty range of manually generated outlines for debris-covered glaciers.

7. CONCLUSIONS AND OUTLOOK

This study has proposed an automatic processing chain to delineate debris-covered glacier boundaries. The processing using InSAR coherence, slope and morphological operations is able to overcome the manual input in conventional optical and thermal remote sensing. In addition to L-band data used in previous studies, we tested the processing chain using X-, C- and L-band radar frequencies along with different processing parameters (spatial resolution, coherence window size, topographic phase). We have shown that C-band Sentinel-1 data without multi-looking and topographic phase removal proved to be the best dataset for our study site. Nevertheless, it should be considered that better results could be achieved at X-band with shorter temporal baselines. We recommend InSAR coherence as a basis for further investigations over stagnant and flowing debris-covered glaciers as data are able to provide valuable support to global databases. Future work should apply our automatic processing chain to map outlines of glaciers located in different regions and climatic conditions. The almost global coverage by Sentinel-1 data, its free accessibility as well as the acquisition of dense time series, make this sensor a suitable base for studies also with shorter temporal baselines and give the possibility for regular updates of global glacier databases in an automatic way. However, more detailed investigation in regard to the acquisition time and temporal base line are recommended for different study regions.