1.

Introduction

The variations in glaciers, ice caps, and ice sheets are considered to be among the sensitive indicators of local or global climate change.^1–4 Most glaciers in the world have shown a general tendency to retreat since the end of the last century.^5–7 Monitoring their evolution is a key issue as the melting of all the glaciers and ice sheets in Antarctica, Greenland, and the high mountains of the world would significantly contribute to the ongoing sea-level rise.^8–11 It is generally agreed that the current climate change poses a threat to the existence of glaciers, especially in low-latitude mountainous regions. If the present warming trend continues, the glacier-covered area will decrease, and most of the glaciers will become thinner or even disappear.^10,12 Glacier variations provide important information associated with global climate change and are also quite helpful in understanding the characteristics of glacier dynamics.^13–15

The displacement field of the glacier surface is the main consequence of glacier activity, and it is essential to fully understand this surface in order to understand the dynamics of an alpine glacier. Such displacement measurements contribute to better knowledge of the rheological parameters controlling the flow of glaciers. These measurements are important in monitoring icefalls, glacier surges, and glacier hazards. They can also detect ice-velocity changes associated with climate change,¹⁶ but it remains difficult and time-consuming to perform regular ground-based surveys of glacier flow even with the advent of differential global positioning system. The field investigations of glaciers provide the most detailed, uncompromised, and reliable information on glacier mass balance and on variations in the glacier size and front, but alpine glaciers are mostly located at high altitudes in remote areas, which significantly increases the difficulties in field surveys. Furthermore, adverse weather, which frequently occurs in such regions, will present additional difficulties in the field. Glaciated areas, therefore, present harsh conditions for field measurements, which means that only a few glaciers have been well documented by fieldwork and even fewer have been monitored with time series of data longer than 10 years.^17,18 Moreover, most of the glaciers that have been monitored are relatively small, simple glaciers that do not fully represent the evolution in local ice conditions, never mind that of the whole world. Also, glaciers and ice in different regions vary in their evolution characteristics.⁵ Because glaciers are numerous, widespread, and mostly remote, remote sensing techniques are needed to acquire comprehensive, uniform, and frequent global observations of glacier variations.^19,20

The observation of glacier changes is coordinated by several international programs, including the U.S. Geological Survey (USGS)-led project, Global Land Ice Measurement from Space, which aims to produce a global glacier inventory from advanced spaceborne sensors. This inventory will enable the study of big questions related to future glacier hydrology and global climate.^19,21 A large amount of archived data has been collected since the launch of the first remote sensing satellite, and this provides a valuable resource for monitoring the long-term changes of glaciers and ice sheets. Satellite remote sensing data were first used to measure glacier and ice sheet changes in Greenland and Antarctic in the late 1970s. Since then, more and more satellite images have been used in systematic large-scale measurements of the state and extent of glaciers and ice sheets all over the world. Older Landsat images, declassified images, aerial-photos, and maps extend our knowledge of glacier changes back to cover most of the last century.¹⁹

Supra-glacial debris exhibits the same spectral properties as the bedrock that lies outside the glacier margin and is thus not detectable by means of multispectral classification alone.²¹ Although supra-glacial debris causes difficulties in automatic glacier delineation, it also produces features on the glacier surface and gives it some of its characteristics—this is helpful for monitoring glacier surface motion using remote sensing data. A detailed map of the ice-velocity field on a mountain glacier can be obtained by cross-correlating optical images, which offers an alternative and effective way of extracting glacier velocity. It can track the displacement of features such as crevasses or surficial debris moving with the ice. It does not require ground control points because the displacement can be referenced to the stationary area. Although the mass balance parameter is difficult, sometimes impossible, to obtain using remote sensing,¹⁹ exploiting glacier surface features in this way can still yield powerful measurement tools for studying glacier state and movements.

The performance of the optical approach has been compared with that of synthetic aperture radar (SAR) interferometry in mountainous regions.²² SAR is another method often used in the study of glacier dynamics.^23–25 Basic interferometric synthetic aperture radar (InSAR) only measures the projection of the displacement vector on to the radar satellite line of sight. SAR data sometimes cannot successfully detect significant displacement of the glacier because of temporal decorrelation.^26,27 Furthermore, the unwrapping operation needed by InSAR also introduces additional noise and decreases the accuracy of the final result. Recently, intensity-tracking and coherence-tracking—two cross-correlation techniques—have been applied to SAR data.²⁸ However, in mountainous regions with significantly rough terrain, SAR images suffer from heavy geometric distortion, such as heavy layover, shadow, and foreshortening, which severely limit the use of this method. Although optical sensors cannot see through cloud and depend on illumination to work, there are a number of operating optical satellites and also a vast amount of archived data that can make up for this deficiency and provide a variety of options for specific applications. The performance of the feature-tracking method has been proofed in applications, and its accuracy and disadvantages compared with both optical and SAR data have been investigated.²⁹ In contrast to the InSAR technique, correlation of optical images directly provides the two horizontal components of the displacement vector. Furthermore, the measurement is unambiguous: absolute displacement can be referenced to stationary areas, which are always present in alpine glaciers. The development of automatic feature-tracking algorithms has significantly increased the accuracy and efficiency of the glacier monitoring approach.^30–32 All these elements make it possible to obtain accurate and detailed measurements of glacier surface motion with Landsat images.

The aim of this study is to demonstrate that a series of surface displacement maps can be routinely obtained on mountain debris-covered glaciers using optical images (Landsat–5/7 band 3 data), and that this approach is more useful than the outline delineation method for monitoring alpine glacier activity. The paper also presents the results of monitoring the movement of the Kuksai Glacier over the past decade, based on fully exploiting the numerous archived Landsat images. This may provide another excellent way of studying glacier evolution in response to climate change.

2.

Study Area and Data Sets

2.1.

Kuksai Glacier

The Kuksai Glacier is located in the Mt. Muztagh Ata region, western China (Fig. 1). Mt. Muztagh Ata [7546 m above sea level (a.s.l.)] is a rounded fault-block mountain at the southeastern edge of the Pamir Plateau (38°15′N, 75°8′E). Several valley glaciers originate near the center of the mountain, flowing away from a central ice cap around the peak. The area is dominated by westerly winds and has a cold, dry climate. According to the meteorological records of Taxkorgan station (37°47′N, 75°14′E, 3090.9 m a.s.l), the average annual and summer temperatures for the period 1957 to 2010 were 3.4 and 15.1°C, respectively, and the average annual precipitation was 70.2 mm.

Fig. 1

(a) Landsat image (bands 743 combined) of Mt. Muztagh Ata, with superimposed elevation contours. The location of this area within China is shown in the inset. (b) Photograph taken from one of the debris-covered glaciers.

Most glaciers in this area are covered by debris to varying degrees. The glaciers are covered by permanent snow and ice above 4700 m a.s.l; below 4700 m a.s.l they are mainly covered by debris. The debris thickness extends from several centimeters at 4700 m a.s.l to several meters at the terminus of the glacier. Generally, the thickness increases with decreasing altitude, even though the debris distribution is affected by several factors such as topography, precipitation, climate, etc.³³ The particle sizes in and distribution of any supra-glacial debris are highly variable, ranging from a few millimeters to blocks several meters across and from minimal occurrence to complete coverage of the glacier tongue.¹⁹ Very small particles are usually distributed uniformly in the ablation area, and the larger-size fractions are often arranged to form characteristic medial moraines. Both of these types of debris would result in a decrease in the reflectivity of the glacier surface, providing appropriate conditions for motion monitoring of the glacier surface with optical data. The supra-glacial debris, most of which has a thickness $>2\mathrm{cm}$, mainly acts as surface insulation for the covered glacier and protects the ice from melting at low altitudes.^34,35 Thus, the underlying surface will contain little glacial melt water, making the glacier surface more resistant.

The ablation season of the Muztagh glaciated region generally starts in May and ends in September. The most rapid snow and ice melting occurs during the months of June to August and reaches its maximum in July.³⁶ The mass balance of the Muztagh glaciers is very poorly sampled in the field and only few glaciers are under regular monitoring.³⁷ In addition, runoff generated by the melting of these glaciers on the mountain is an important source of water for people living downstream and also for hydropower stations. So information on the glaciers is also necessary when planning the development of towns and the construction of hydropower stations. The ongoing glacier retreat, therefore, has important social and economic impacts, and measuring this retreat is a first step toward the prediction of future water resources in the glaciated area.¹⁰

Covering $>86.5{\mathrm{km}}^2$, the Kuksai Glacier constitutes an important component of the Muztagh glacier area. The Kuksai Glacier is located in a large valley on the east side of Mt. Muztagh, with an elevation range of 4000 to $\sim 5600\mathrm{m}$ a.s.l. Its total length is $\sim 18\mathrm{km}$ and has a maximum width of 1.5 km. Because of the relatively low gradient, the weight of the ice provides only a small force to make the ice flow downward. There are two sub-glaciers in the upper part, which is the main area of accumulation, and these converge into one at 4600 m a.s.l (see Fig. 1). Given the size and remoteness of the Kuksai Glacier, satellite imagery is a suitable means of obtaining a comprehensive and more regular sampling of it.

Remote sensing studies based on advanced spaceborne thermal emission and reflection radiometer images and old aerial photographs have concentrated on observing the extent of glaciers in the Mt. Muztagh Ata area.³⁸ The results show that most of the glaciers have been retreating during the observation period. Because of the increased melting of snow and ice with the warming temperatures, the meltwater production has increased and often accumulates in closed basins on the glacier surface. Therefore, some moraine-dammed lakes on the tongue of the glacier have been detected from remote sensing images. Although the lake area is small, the appearance of the supra-glacial lake is significant proof of changes in the local climate.³⁹

3.

Methods

3.2.

Computing the Motion of the Glacier

The motion of a glacier can be computed by tracking the displacement of features in two optical images of the same area taken at different times, based on the NCC method. For some outlet glaciers of the Greenland or Antarctic ice sheets, the persistence of the surface features permits velocity measurements from images separated by as long as 11 years. Some velocity fields have also been derived from optical images separated by more than a year on mountain glaciers.²⁷ The low contrast on the glacier surface is an obstacle to deriving the flow map. In order to overcome the inaccuracy associated with the low contrast, the preprocessing step of image intensity enhancing is introduced into the data processing.

The key step in computing glacier surface motion is a precise relative orientation of two images obtained by adjusting the slave image, assuming that the master image is well georeferenced. This is performed using numerous precisely matched points extracted automatically based on the NCC method.²⁷ The NCC algorithm is an intensity-based similarity measure, which is used in image matching to measure the similarity between matching entities in one image and their corresponding entities in another image. The algorithm was developed based on the concept of distance measurement but second normalized to account for the differences in brightness and contrast.⁴¹ The NCC coefficient is computed according to Eq. (1).

Eq. (1)

$$\mathrm{NCC}(u,v)=\frac{\sum _{x,y}[f(x,y)-{\overline{f}}_{u,v}][t(x-u,y-v)-\overline{t}]}{{\{\sum _{x,y}{[f(x,y)-{\overline{f}}_{u,v}]}^2\sum _{x,y}{[t(x-u,y-v)-\overline{t}]}^2\}}^{0.5}},$$

where $\mathrm{NCC}(u,v)$ $a$ is the correlation coefficient, $f$ is the search window in the slave image, and $t$ is the counterpart template window in the master image. ${\overline{f}}_{u,v}$ is the mean of the subwindow $t^{\prime }$ of the search window $f(x,y)$ underneath the template whose top-left corner lies at pixel ($u,v$), as shown in Fig. 2. $\overline{t}$ is the mean of the template windows.

Fig. 2

The normalized cross-correlation scheme.

The search window, $f$, is usually bigger than the template window, $t$. $\mu$ and $\nu$ are the offsets between the two windows in the $x$- and $y$-direction, respectively. The Euclidean distance ($\mu ,\nu$) between the coordinates of the reference point $[x,y]$ and the matching point $[x-u,y-v]$ is considered as the displacement. Even if there is no truly corresponding entity in the search window, there will always be some peak correlation coefficient, which should be discarded. It should also be noted that a small NCC coefficient could be related to a feature variation in the selected areas during the time period considered.³¹ Therefore, it is necessary to decide on a threshold of $\gamma$ for the NCC coefficient, below which the match should be rejected. In this study, we empirically set the threshold value equal to 0.3, which gave satisfactory results.

Only the pixel-level accuracy peak locations can be obtained directly from the original, uninterpolated images based on the NCC coefficient. However, those locations might not be the exact positions of the matching entities. In order to obtain the positions with subpixel accuracy from the discrete NCC coefficients, in this paper, coregistration subimage pair (template and search window) interpolation based on the fast Fourier transform method is performed with an oversample factor of 100. In addition, it is necessary to match images with the same track and frame number because the different sun elevations will give a variation shadow, which presents an obscuring factor for the image processing. The size of the template should be chosen to be large enough to maximize the signal-to-noise ratio and small enough to minimize velocity gradients.⁴² The search window should be chosen to be large enough to include the farthest moving template and small enough to limit the computational cost of the matching.³¹ Based on empirical considerations, a template widow size of $32\times 32\text{pixels}$ and a search window of $64\times 64\text{pixels}$ was used in this study, and the original coregistered Landsat TM image pair was first oversampled twice. The sample window size was set to be half that of the template window. These parameters can keep the resulting flow field smooth and also preserve the local features of glacier motion. The surface velocity field measurement method based on NCC offers an improved method of retrieving information about the behavior of the studied glacier. The glacier deformation varies with the seasons and the local slope, but most features on the glacier surface remain stable for a relatively long time, which is the key prerequisite for using the NCC method. Potentially, this method can also be used to derive displacements of the Earth’s surface, especially when based on higher spatial resolution images.⁴³

4.

Results and Analysis

Several Kuksai Glacier motion maps, see Fig. 3, were derived from the multitemporal Landsat image data. The general trend of the glacier surface displacement is similar on the whole: the magnitude of the glacier motion is high in the upper part of the glacier and decreases with decreasing altitude. According to the differences in the motion distribution pattern, the debris-covered Kuksai Glacier can be generally divided into several subregions, such as the downstream, middle, and upstream parts. There are also some differences between the velocity maps for different years in the same subpart.

Fig. 3

The flow maps for 1998 to 2003 and 2009 to 2010 (unit: meter) superimposed on the Landsat TM band 3 data. (a) Flow components of the glacier along the north–south and east–west directions. (b) The overall glacier flow.

The glacier flow maps show an apparently similar movement pattern over a period of approximately one year, generally including one ablation season and one accumulation season. Statistically, the NCC describes the correlation between images acquired at different times, and so it can be considered as a function of the time difference. It is possible to hypothesize that a mountain glacier flows at a low velocity during the period being considered. In order to avoid the signal of the glacier motion being buried by noise in the results, it is important to select images with sufficient time separation. One year is enough for efficiently extracting the glacier motion signal. It is also convenient for comparing the results between different years. The ice stream behaves differently in the different portions of the glacier, and the different velocities are shown in different colors in Figs. 3(a) and 3(b).

The maps in Fig. 3(a) show the displacement of the Kuksai Glacier in the north–south (NS) and east–west (EW) directions in the years 1998 to 2003 and 2009 to 2010. Some obvious features can be seen in the NS maps; these are due to the shape and distribution of the glacier. The motion of the upstream and downstream parts of the glacier apparently has large components along the NS direction but with opposite orientations. In the NS maps, the middle section appears to be motionless because the direction of motion and the NS axis are perpendicular to one another. However, the corresponding part of the glacier, colored dark red in the EW maps in Fig. 3(a), has a high degree of movement along the EW direction. The motion of the glacier surface, see Fig. 3(b), gives an excellent description of the motion characteristics of the Kuksai Glacier. The average approximate annual glacier motion is clearly presented, although it is influenced by noise. The map for 2000 to 2001 suffers from more noise than other years in the middle and upstream parts because of the poor quality of the Landsat TM data. This resulted from atmospheric contamination due to cloud cover in the upper part of the Kuksai Glacier. Due to their covering a much longer part of the ablation season (August 16, 1999, to September 11, 2000), the NS and EW maps show a much larger deformation, especially in the EW direction, which is in great contrast with the result for 2002 to 2003. Because of a lack of appropriate Landsat image data, the flow map of the 2002 to 2003 period shows much less motion than the others owing to the short time interval, which is only about eight months. Furthermore, it covers very little of the main part of the ablation season of 2002 to 2003, and so only a small movement can be observed in the glacial area, as seen in Figs. 3 and 4. The southern branch of the upper part of the glacier has few useful match points because it was always covered by snow and because of a lack of available features in the corresponding part of the image data. So, this part is not covered by the profiles shown in Fig. 4. The polynomial approximation of motion (see Fig. 4) in the center streamline was computed based on all the data except that for the period 2002 to 2003 as this cannot fully represent the annual motion characteristics of the Kuksai Glacier for that year.

Fig. 4

Variation in glacier motion at the center streamline as a function of distance from the front and polynomial line fit (excluding data for 2002 to 2003).

From the motion distribution maps, it can be clearly observed that the ice flow vector varied along with position on the surface of the glacier. The spatial variation is an obvious reflection of the glacier variable activities and provides information on the state of the glacier, which can be used in detailed studies of the glacier dynamics. At the glacier front, as shown in Fig. 3(b), the ice appears to be stable and has no apparent motion. This result is highly consistent with previous studies that used the outline delineation method. From Fig. 3 it can also be found that the flow of the glacier’s edges is slow relative to that of the center, which is because of the friction caused by the valley’s side walls. A similar phenomenon has also been observed in other studies.^27,44,45

The motion is approximately constant from 9 to 11 km upstream of the glacier front, from where the glacier flow begins to decelerate until the front. The fluctuation of the glacier motion between different years in the middle and upper parts is much larger than that in the downstream part, which can be clearly observed in Fig. 4. There is a special location 11 to 13 km upstream from the front, which is the convergence region of the two upstream branches. Qualitatively, we observed a conspicuous acceleration upstream of the intersection and also a slight but distinct acceleration downstream of the intersection in the fitting curve.

The maximum magnitude of the glacier flow was found to be $50\mathrm{m}/\mathrm{a}$ in the upper part of the Kuksai Glacier. There are two portions with relatively high velocities in the study glacier: the first is the middle part (about $45\mathrm{m}/\mathrm{a}$) and the other is the northern branch of the upstream part ($50\mathrm{m}/\mathrm{a}$). However, glacier displacement maps lack information on the origin of the glacier because of the snow cover and lack of matching features. It can be seen that the velocity of the downstream section is roughly constant over the section 0 to 2 km from the front. This situation has continued largely unchanged over the 12 years since 1998, judging from the profiles obtained from the flow magnitude map derived from Landsat observations shown in Fig. 4.

The velocity distribution pattern of the ice flow shows a few variations in different years. It can be seen that the motion varies considerably from year to year 8 to 16 km upstream from the front. There is a marked average difference in the flow magnitude of $\sim 10\mathrm{m}/\text{year}$ between 1999 and 2000, and 2009 and 2010. Compared with the middle and upstream parts, the average variation in the glacier motion between different years in the downstream part of the glacier is small. The approximate annual velocity in 1998 to 2003 was larger than that in the period 2009 to 2010 under the condition of equal time intervals. The approximate annual velocity became a little slower with time, see Fig. 4. The velocity fluctuations in the upstream part are much larger than those in the downstream part, which means that the most intense activity occurred in the upstream part of Kuksai Glacier.

5.

Flow Velocity Error Analysis

There are many surface features captured in image scenes, such as mountain rock, glacier, road, artificial construction, and so on. They can be divided into two categories, active and inactive, based on their motion. In our study area, the mountain rocks around the glacier and the most part of the remote sensing image except the glacier-covered region are to be considered as stable, which provide a valuable reference for determining the glacier velocity.

As motioned previously in Sec. 3, the template/search window and the image spatial resolution will heavily affect the accuracy of the last result. The NCC method used in this study can give the theoretical registration error $<0.1\text{pixel}$. But due to the presence of error in data processing, the reduction of image quality caused by bad weather, and the contribution of the residual topography, altitude, and satellite orbit, the final error will be $>0.1\text{pixel}$. There is another important error factor originating from the scan operation mode of satellite. The Landsat satellites adopt the whisk scan operation mode to scan the earth surface and put the many scan strips together for scene data. Due to the data acquired mode, there will be a little alignment error between two adjacent scans. So this error will also contribute to the total error of the final result.

Displacement can be referenced to motionless areas, which are always available for mountain glaciers.²⁷ The noise from both the scan operation mode and the registration method is directly reflected as the residual value in the nonmotion area. In order to obtain the value of the error, we can statistically analyze the residual error in the non-glacier region, especially in the area without rough topography plane. The accuracy of results is 7 m with corresponding standard deviation of 5 m. The error caused by the sensor operation mode is also included. The way of error estimation in the method is based on the motionless area, which usually is non-glacier region.

6.

Conclusion

Retrieving the fluctuations and movement of the Kuksai Glacier, the biggest glacier in the Mt. Muztagh Ata area, using archived Landsat imagery was the main objective of this paper. In this study, the glacier surface displacement information was efficiently retrieved from the image sequences of the archived optical Landsat satellite imagery. The flow vectors of the Kuksai Glacier surface were found using image pairs in which the images were taken nearly 1 year apart and using the NCC method based on images obtained over a total span of 6 years (1998 to 2003 and 2009 to 2010). The flow measurement provides relatively detailed glacier surface flow information, which is an extremely valuable parameter in glacier dynamic studies. The ice displacement maps play an important role in understanding the evolution of the glacier and its response to climate change. There is no apparent activity in the Kuksai Glacier front in all the extracted velocity results. Its terminus has proved to be stationary over the past decade.

The general trends for the Kuksai Glacier are similar to those of some other large mountain glaciers: the ice flows faster in the upstream part than in the downstream part. There is also no big difference in the maximum velocity (usually less than $\sim 20\mathrm{m}$ per year) or in the general flow patterns between the Kuksai Glacier and other large alpine glaciers in the Himalaya and Alaska,^44,46,47 which also effectively reflects the motion characteristics of large mountain glaciers.

The detailed study of year-by-year variations in glacier behavior using the NCC method still can reveal a significant change in the velocity of the upstream portion of Kuksai Glacier between different years. It is obvious that the fluctuation in the velocity in the upstream and middle parts is larger than that in the downstream part.

7.

Discussion

The gentle terrain in the front of Kuksai Glacier and the thick debris-cover cause the ice in the front to cease to move and to melt in situ. This presents appropriate conditions for the appearance of dead ice, which is usually distributed in the glacier terminus. The information regarding the motion of the glacier front in the Kuksai Glacier distinctly reflects the existence of dead ice. It heavily restrains the terminus retreat of the Kuksai Glacier and keeps the glacier terminus relatively stable. So there is no apparent movement or variations in glacier outline at the glacier front over a long period. The low correlation between glacier extent and climate reinforces the need for measurements of glacier movements. In order to study the detailed state of the Kuksai Glacier, it is essential to estimate the glacier surface flow field.

There are few available features above the snow line because of the snow cover on the glacier surface all year round. Below the snow line, the glacier surface is covered by debris, which gives more features for extracting motion information based on the optical image cross-correlation method. Therefore, the obtained glacier flow fields also usually only cover the part of the glacier surface below the snow line. It can be observed that the ice velocity generally becomes smaller with decreasing altitude in all the obtained flow fields below the snow line.

The confluence of the big southern sub-glacier makes the motion pattern of the Kuksai Glacier complicated. The sub-glacier provides an important ice source for the Kuksai Glacier. The supplement from the sub-glacier increases the ice quantity, which stops the decrease in the flow velocity and, in fact, causes it to increase. Then, after reaching a maximum, the velocity decreases again. The increase of the ice resource leads to the corresponding velocity variation because of the new mass supplement from the sub-glacier.

From six flow field maps of the glacier, a detailed analysis of glacier surface motion trajectory patterns may allow for the detection and measurement of tidal effects in the glacier motion, which is another advantage of the method described over the traditional method. The fluctuations in the glacier movement between different years in the middle and upper parts, especially in the section 9 to 16 km upstream from the glacier front, are much larger than that in the downstream part, which may be related to the variation of glacial thickness and climate elements such as precipitation, temperature, and wind. Further research is needed to illuminate the details and reasons for this behavior.