1. INTRODUCTION
The accurate knowledge of glacier extent is a prerequisite to many glaciological studies,for example volume change estimates (Kääb and others, Reference Kääb, Berthier, Nuth, Gardelle and Arnaud2012; Gardelle and others, Reference Gardelle, Berthier, Arnaud and Kääb2013) and ice dynamic modeling (Huss and Farinotti, Reference Huss and Farinotti2012; Fürst and others, Reference Fürst2016). Change rates of glacier front positions are often used variables in studies of climate change (Ranzi and others, Reference Ranzi, Grossi, Iacovelli and Taschner2004; GCOS, Reference GCOS2011) or natural hazards such as glacial lake outburst floods (Bajracharya and Mool, Reference Bajracharya and Mool2009; Hewitt, Reference Hewitt2014). 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, Reference Ranzi, Grossi, Iacovelli and Taschner2004; Pfeffer and others, Reference Pfeffer2014; Paul and others, Reference Paul2015). 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, Reference Paul2015; Smith and others, Reference Smith, Bookhagen and Cannon2015) and thermal (Alifu and others, Reference Alifu, Tateishi and Johnson2015) 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, Reference Alifu, Tateishi and Johnson2015). Optical data are often limited due to the lack of cloud-free acquisitions and good illumination (Smith and others, Reference Smith, Bookhagen and Cannon2015). 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, Reference Paul, Huggel and Kääb2004). However, Smith and others (Reference Smith, Bookhagen and Cannon2015) improved the algorithm of Paul and others (Reference Paul, Huggel and Kääb2004) 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, Reference Bolch, Buchroithner, Kunert and Gomarasca2007; Shukla and others, Reference Shukla, Arora and Gupta2010; Bhambri and others, Reference Bhambri, Bolch and Chaujar2011; Racoviteanu and Williams, Reference Racoviteanu and Williams2012). 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, Reference Schauwecker2015). 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, Reference Rodriguez and Martin1992; Michel and Rignot, Reference Michel and Rignot1999; Li and others, Reference Li, Guo, Li and Wang2001; López-Martinez and Pottíer, Reference López-Martínez and Pottier2007) 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, Reference Atwood, Meyer and Arendt2010), the Pamirs (Jiang and others, Reference Jiang, Liu, Wang, Lin and Long2011) and the Himalaya (Strozzi and others, Reference Strozzi, Paul and Kääb2010; Frey and others, Reference Frey, Paul and Strozzi2012). Atwood and others (Reference Atwood, Meyer and Arendt2010) coupled the coherence with slope, size and morphological filters to map glacier outlines. However, coherence images of Himalayan glaciers, used in Frey and others (Reference Frey, Paul and Strozzi2012), contained data voids due to steep terrain. Therefore, additional manual editing was required. Wu and others (Reference Wu, Zhang, Zhang, Lu and Zhong2012) 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 (Reference Robson2015) combined InSAR coherence with object-based image classification (Rastner and others, Reference Rastner, Bolch, Notarnicola and Paul2014) 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, Reference Frey, Paul and Strozzi2012; Smith and others, Reference Smith, Bookhagen and Cannon2015). 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 km2 (Rankl and others, Reference Rankl, Kienholz and Braun2014). 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, Reference Bajracharya, Shrestha, Bajracharya, Maharjan and Guo2014). 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, Reference Bajracharya, Shrestha, Bajracharya, Maharjan and Guo2014) and Rankl and others (Reference Rankl, Kienholz and Braun2014) only partially integrate debris-covered parts into their glacier outlines, whereas the RGI 5.0 glacier outlines (Cogley and others, Reference Cogley2015) 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 (Reference Rankl, Kienholz and Braun2014).
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, Reference Weydahl2001). 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, Reference Atwood, Meyer and Arendt2010; Weydahl, Reference Weydahl2001). As reference data, we used outlines from the updated global glacier inventories GLIMS (Bajracharya and others, Reference Bajracharya, Shrestha, Bajracharya, Maharjan and Guo2014), the RGI version 5.0 (Cogley and others, Reference Cogley2015) and the glacier outlines from Rankl and others (Reference Rankl, Kienholz and Braun2014) (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 (Reference Rankl, Kienholz and Braun2014) 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, Reference Michel and Rignot1999; Frey and others, Reference Frey, Paul and Strozzi2012). The coherence of two images S 1 and S 2 can be calculated by averaging over M × N pixels (coherence window) in the image dimensions m × n, mathematically defined as
$$\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, Reference López-Martínez and Pottier2007). 
$S_{2}^{\ast}$ is the complex conjugate of S 2.
If other decorrelation effects like thermal noise, geometric or volume decorrelation (Li and others, Reference Li, Guo, Li and Wang2001; Atwood and others, Reference Atwood, Meyer and Arendt2010) 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, Reference Werner, Wegmüller, Strozzi and Wiesmann2005). 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, Reference GAMMA2011). 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, Reference Veci2015). 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, Reference Touzi, Lopes, Bruniquel and Vachon1999). 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, Reference Frey, Paul and Strozzi2012). 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, Reference Lee and Liu2001). Therefore, previous authors suggested subtracting the topographic phase from the complex interferogram in order to account for such effects (Atwood and others, Reference Atwood, Meyer and Arendt2010; Jiang and others, Reference Jiang, Liu, Wang, Lin and Long2011). 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 (Reference Atwood, Meyer and Arendt2010), coherence and slope thresholds were applied. Bolch and others (Reference Bolch, Buchroithner, Kunert and Gomarasca2007) reported a maximum average slope of 10° for the debris-covered parts of the Mount Everest region in ASTER DEM. Paul and others (Reference Paul, Huggel and Kääb2004) 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 (Reference Atwood, Meyer and Arendt2010) 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, Reference Haralick and Shapiro1992). 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, Reference Strozzi, Luckman, Murray, Wegmüller and Werner2002; Seehaus and others, Reference Seehaus, Marinsek, Helm, Skvarca and Braun2015). 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 (Reference Burgess, Forster, Larsen and Braun2012). 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 (Reference Wegmüller and Werner1995) 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 (Reference Huang, Tian, Li and Zhou2017) 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 km2 in Fig. 4f and 0.88 km2 in Fig. 4d) than Sentinel-1 (1.69 km2, 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 km2. 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−1, whereas the upstream glacier flow reaches up to 0.60 m d−1 (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, Reference Holzer2015; Vijay and Braun, Reference Vijay and Braun2016; Ragettli and others, Reference Ragettli, Bolch and Pellicciotti2016). 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, Reference Santoro, Werner, Wegmüller and Cartus2007) 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 (Reference Alifu, Tateishi and Johnson2015), 2–10% in Smith and others (Reference Smith, Bookhagen and Cannon2015), 11.5% in Rastner and others (Reference Rastner, Bolch, Notarnicola and Paul2014), 7–17% in Robson and others (Reference Robson2015)). 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, Reference Frey, Paul and Strozzi2012). 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.
ACKNOWLEDGMENTS
This work was financially supported by the project BR 2105/14-1 within the DFG Priority Program ‘Regional Sea Level Change and Society’, the DLR/BMWi project TanDEM-ICE (FKZ 50E1414) as well as the HGF Alliance ‘Remote Sensing and Earth System Dynamics’ (HA-310). TerraSAR-X data were ordered from DLR through AO HYD1788 and LAN0164; ALOS data via the ESA AO 3575 and via ASF. We thank for cost-free distribution of Sentinel-1 data via the Sentinel Open Access Hub as well as USGS for free access to Landsat imagery via EarthExplorer. Special thanks go to Dr. Jenny Turton for her support as a native speaker in the improvement of the manuscript. We acknowledge support by Deutsche Forschungsgemeinschaft and Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU) within the funding pogramme Open Access Publishing.
Open access
