Accurate Determination of Glacier Surface Velocity Fields with a DEM-Assisted Pixel-Tracking Technique from SAR Imagery

Abstract

We obtained accurate, detailed motion distribution of glaciers in Central Asia by applying digital elevation model (DEM) assisted pixel-tracking method to L-band synthetic aperture radar imagery. The paper firstly introduces and analyzes each component of the offset field briefly, and then describes the method used to efficiently and precisely compensate the topography-related offset caused by the large spatial baseline and rugged terrain with the help of DEM. The results indicate that the rugged topography not only forms the complex shapes of glaciers, but also affects the glacier velocity estimation, especially with large spatial baseline. The maximum velocity, 0.85 m∙d⁻¹, was observed in the middle part on the Fedchenko Glacier, which is the world’s longest mountain glacier. The motion fluctuation on its main trunk is apparently influenced by mass flowing in from tributaries, as well as angles between tributaries and the main stream. The approach presented in this paper was proved to be highly appropriate for monitoring glacier motion and will provide valuable sensitive indicators of current and future climate change for environmental analysis.

1. Introduction

Glaciers, especially mountain glaciers, are usually treated as one of the most sensitive indicators of the Earth’s climate change, which has always been characterized by natural variability such as global annual air temperature and average precipitation [1,2,3,4]. The former shows significant increase during the last century [5,6]. Nowadays, the great majority of glaciers experience on-going retreat, particularly mountain glaciers in the mid- and low-latitudes, such as the Himalayas, Alps, and Pamirs [1,2,7,8]. This phenomenon was more serious at the end of last century and was identified in most of the examined glaciers on the Tibetan Plateau and in Europe, although a small number of glaciers continue to advance [1,9,10,11]. In addition, ice melting in response to global warming in glaciated areas could further accelerate glacier flow and potentially lead to significant loss of glacier mass, which would affect river runoff, cause glacial retreat or bring about supraglacial lake outburst floods [7,9].

Extracting the surface movement of mountain glaciers plays a significant role in investigating the glacier mass balance, because the ice mass flow is highly dependent on the glacier flow velocity [12,13,14]. Therefore, knowledge of the glaciers’ surface velocity fields is highly significant in understanding the response of glaciers and ice caps to climate changes and predicting glacier-related hazards [3,15]. Given the poor spatial sampling and difficulty of monitoring glaciers directly in remote regions by fieldwork, satellite remote sensing presents an attractive alternative for efficiently monitoring velocities with optical and radar imagery [9,16,17,18]. However, optical remote sensing for glacier monitoring is significantly limited by the unfavorable weather conditions of glacial areas, such as thunderstorms, rain and snowfall, all of which usually cause frequent and rapid changes in illumination. Additionally, it is difficult to routinely obtain cloudless optical data for these regions, which is required for regular monitoring of glacier flow distributions.

Compared to optical remote sensing methods, Synthetic Aperture Radar (SAR) has the ability to image through cloud cover and work day-and-night. It is insensitive to variations in weather or illumination. It can obtain accurate results efficiently by extracting the phase difference of SAR imagery. Therefore, the interferometric SAR (InSAR) method is a valuable technique for studying glacier dynamics due to its high sensitivity to ground deformations [19,20,21]. Unfortunately, the InSAR method frequently suffers from the heavy decorrelation caused by scatter property changes over a period or the steep gradient of motion associated with fast ice movement [22,23].

However, the pixel-tracking approach based on SAR intensity information can still be used for extraction of ice motion. The excellent advantage of pixel-tracking over InSAR is the fact that pixel-tracking is hardly influenced by temporal and spatial coherence. The strict limitations on the time interval and spatial baseline of InSAR are no long necessary, as long as the SAR imagery is from the same area [16,24]. Furthermore, the pixel-tracking method can yield two-dimensional motion from a SAR data pair, rather than only the movement along the line of sight (LOS) of the SAR sensor. Therefore, it is an excellent alternative to the InSAR technique, especially in high-accumulation or high-melt glaciers in mountainous areas [4,22,25]. Based on C-, L- and X-band SAR imagery, it has been successfully applied to calculate the velocities of glaciers in several diverse locations [4,15,22,26,27,28], e.g., Greenland, Antarctica, the Arctic, and mountain glaciers in alpine regions. But, the complex, rugged terrain of mountainous glacial areas usually yields an external component not related to the glacier flow signal, which should be estimated and removed for accurate ice motion.

In this study, we used the refined pixel-tracking method to efficiently derive glacier surface velocity fields with SAR intensity images acquired over glaciers in Central Asia using ALOS/PALSAR imagery, which has a nickname Asia’s water tower because of providing water storage in the large proportion of ice cover in this region. The objective of this work was to calculate the accurate and detailed velocity distribution of these mountain glaciers efficiently and analyze their spatial flow characteristics.

2. Study Area and Datasets

The Pamir Mountains, with their dense glacier distribution, are mainly located in the Central Asian. The study area was mostly in the southern portion of Tajikistan, covering the geographic area between 38°15′–39°20′N and 72°00′–72°40′E. The altitude generally varies from 2800 m to 7400 m above sea level (a.s.l.). The highest elevation in this region is the Ismoil Somoni Peak, with an altitude of 7495 m a.s.l. There are many other famous high mountains around the region, including Garmo Peak and Independence Peak. All of these contribute to the complex terrain and deep valleys in this mountainous glacial area. Consequently, it provides the proper generation and development conditions for large mountain glaciers in the Pamirs.

Unlike precipitation in the Himalayas and Tian Shan [29,30], precipitation in the Pamirs is usually greatest during the winter and spring due to the Siberian/Tibetan high-pressure systems during the winter [31]. The mean annual temperature and precipitation from the Fedchenko Observatory Station were estimated to be about −6.7 °C and 1200 mm, respectively. The climatic conditions, as well as rugged terrain, provide a perfect precondition for the development of glaciers and also affect their evolution, such as Fedchenko Glacier, Grumm-Grzhimaylo Glacier, Medvezhiy Glacier, Nalivkina Glacier and Bivachny Glacier, whose names were derived from the Global Land Ice Measurements from Space (GLIMS) database and previous studies. We also coded some glaciers in the study area as No. 1–7 because most of them had no published names (Figure 1). Their basic information is presented in Table 1, which is directly measured and obtained from the Landsat-7 imagery.

Table 1. The basic information of some glaciers in study area.

Glacier Name	Length(km)	Area(km2)
Fedchenko	77	146
Grumm-grzhimaylo	34	56
Nalivkina	17	26
Bivachny	28	31
No. 1	7	12
No. 2	11.5	22
No. 3	5.5	7
No. 4	7	11
No. 5	15	23
No. 6	11.5	10
No. 7	12.5	9

Table 1. The basic information of some glaciers in study area.

Glacier Name	Length(km)	Area(km2)
Fedchenko	77	146
Grumm-grzhimaylo	34	56
Nalivkina	17	26
Bivachny	28	31
No. 1	7	12
No. 2	11.5	22
No. 3	5.5	7
No. 4	7	11
No. 5	15	23
No. 6	11.5	10
No. 7	12.5	9

Fedchenko Glacier is about 77 km long and from 1700 to 3100 m wide. It is by far the world’s largest glacier found outside of the polar region and runs along the eastern mountainside of the Science Academy mountain range. The glacier’s elevation ranges from 2900 to 6280 m and its equilibrium line altitude is approximately situated around 4700 m [31]. The Fedchenko Glacier consists of several small tributaries and a valley glacier-type tongue, which extends from south to north in the center of the Figure 1 (light blue on the upstream part and red on the downstream section). Since the upstream glacier covered with snow and ice reflects sunlight quite strongly, it looks light blue, whereas the surface of the downstream section is covered with supraglacial debris, so it looks red as it reflects sunlight rather weakly.

The glacier meltwaters in this region mainly feed into the Surkhob River and Amu Darya River, which supply the agricultural and economic activities downstream of rivers in the mountains of the study area. Many glaciers in the Pamirs, as well as the Fedchenko Glacier, have been experiencing accelerating melting in recent years. Thus, their influence on water resources is gradually becoming a serious concern in these area [31,32,33]. The historical surface velocity observations showed horizontal displacement of 0.5–0.7 m per day in the central part of Fedchenko glacier and decreasing values toward the glacier tongue [33]. The comparatively small fluctuation in ice area proves that glacial outline change detection is inefficient for studying glacier activity, especially within a short temporal separation. But the pixel-tracking method based on SAR intensity information can provide valuable motion information, which is a much better indicator for climate-related glacier reactions.

Figure 1. Landsat-7 image of the study area. The cyan represents the mountain glacier distribution. The yellow outline indicates the coverage of SAR imagery. The white profiles are along the center flowlines glaciers. Capital letters A–F show the junction positions. The area of interest is indicated in the inset map.

Figure 1. Landsat-7 image of the study area. The cyan represents the mountain glacier distribution. The yellow outline indicates the coverage of SAR imagery. The white profiles are along the center flowlines glaciers. Capital letters A–F show the junction positions. The area of interest is indicated in the inset map.

Compared with C- and X-band SAR data, the L-band SAR works at a much longer wavelength (23.6 cm) and provides higher coherence than other SAR data, particularly on dense vegetation and mountain glaciers [34]. L-band can penetrate deeper, maintaining much more robust features on the glacier surface [35]. C- and X-band data often suffer from decorrelation problems due to their shallower penetration depth in snow and ice [15,36], especially if solid precipitation occurred within the period of SAR data acquisitions. Accordingly, L-band SAR data can complement existing applications challenged with C-band and X-band data [25,27].

For obtaining motion distribution on the glacier surface, we needed two SAR images acquired from repeat orbits nearly covering the same area. The stability of the features on the glacier surface is often affected by weather conditions and the movement of ice. Therefore, SAR imagery with 46-day separation was collected from the L-band operational PALSAR sensor loaded on the Advanced Land Observation Satellite (ALOS), which was launched by the Japanese Aerospace Exploration Agency (JAXA) in 2006. The SAR remote sensing imagery used in this study is summarized in Table 2.

Table 2. ALOS/PALSAR image pairs used in pixel-tracking processing.

Date Acquired	Orbit No./Frame No.	Perpendicular Baseline	Temporal Baseline	Polarization	Orbit Type
2 January 2007	530/760	1825 m	46 days	HH	Ascending
17 February 2007	530/760

Table 2. ALOS/PALSAR image pairs used in pixel-tracking processing.

Date Acquired	Orbit No./Frame No.	Perpendicular Baseline	Temporal Baseline	Polarization	Orbit Type
2 January 2007	530/760	1825 m	46 days	HH	Ascending
17 February 2007	530/760

Several high mountain peaks and deep valleys form the complex, rugged terrain in the study area, which may introduce an additional displacement value in the displacement field due to the stereo effect and seriously contaminates the accuracy of the velocity measurement results. For better understanding and precisely eliminating the topographic effect, an external digital elevation model (DEM) of the study area is required. In this study, SRTM DEM version 4 with nominal resolution of 3 arc-seconds or 90 m is employed as it is sufficient for compensating for the topographic contribution.

3. Method

3.1. Pixel-Tracking Approach Analysis

The fundamental hypothesis of the pixel-tracking method for extracting glacier displacement is that the movement only occurs in the glaciated area, whereas the ice-free region is supposed to maintain the stability of quiescence. In other words, the pixel-tracking approach is a relative deformation measurement method like the InSAR technique. Therefore, the SAR intensity based pixel-tracking approach is a promising alternative to SAR interferometry for estimation of glacier motion, because of the serious decorrelation caused by rapid and incoherent flow or changes in the glacier’s surface roughness [22,23]. This technique is mainly based on the normalized cross-correlation (NCC) maximization between SAR intensity image patches from SAR imagery pair respectively, which is a useful and efficient similarity measurement method for automated image coregistration process [16,25]. The patch window with the peak of the correlation coefficient within the search area could be considered to be the most similar window and used to calculate the ice motion [22]. The accuracy of a small fraction of a pixel can be achieved by match peak determination through oversampling of the correlation surface. Additionally, the reliability of the individual motion vector can be determined by the signal-to-noise ratio (SNR) of the correlation coefficient.

Firstly, the SAR Single Look Complex (SLC) data was generated by the Repeat Orbit Interferometry Package (ROI_PAC) software from the raw PALSAR data. Then, the pixel-tracking method was immediately applied to the PALSAR SLC data. The search window should be large enough to ensure that the largest possible displacement was included. A larger patch window will generally yield a higher matching accuracy, although the match procedure will require a longer processing time. Hence, a patch size of 32 × 32 pixels for extracting the match peak with a step of 4 pixels in the range direction and 12 pixels in the azimuth direction was employed in the fine match operation for obtaining an optimal balance between the glacier velocity precision and processing time [24]. Prominent surface features that could be identifiable on two co-registered images, such as crevasses, rifts and glacial edges that move associated with the ice, are also helpful to accurately determine the displacement.

The range offset and azimuth offset given by sub-pixel NCC both consist of the glacier motion-related signal and no-motion component. The latter is mainly contributed by the imaging geometry and rugged terrain in the study area. For precisely obtaining the glacier velocity, the no-motion part in both offset results should be removed. Generally, the components of the offset value retrieved by sub-pixel NCC can be clearly described by the following equation:

$$offse{t}_{total}=offse{t}_{ice}+offse{t}_{orb}+offse{t}_{att}+offse{t}_{ion}+ offse{t}_{ato}+offse{t}_{top}+offse{t}_{res}+offse{t}_{noi}$$

(1)

where $offse{t}_{ice}$ is glacier surface motion related to offset; $offse{t}_{orb}$ is the offset caused by the orbit residual error; $offse{t}_{att}$ is the offset value induced by sensor attitude； $offse{t}_{ion}$ is the offset caused by the ionosphere; $offse{t}_{ato}$ is atmospheric-related offset, which often can be ignored due to its insignificant effect;$offse{t}_{top}$ is the offset value caused by complex terrain; $offse{t}_{res}$ is the offset value induced by global residual error, which is usually caused by errors in global transfer parameter estimation; and $offse{t}_{noi}$ is noise related offset caused by data processing.

Usually, the offset associated with orbit and sensor attitude between two SAR images is a function of range and azimuth position with a second-order polynomial relationship. With this matching-function relation, thus, the components corresponding to the long-wavelength signals including satellite orbit and sensor attitude in the offset fields can be modeled and removed by a second-order polynomial fitting in an ice-free glaciated area. As mentioned above, the atmospheric variation during the interval period of SAR the image pair has little influence to the last displacement field because of its comparatively low accuracy, unlike the InSAR approach, and it can be ignored correctly, especially for the SAR imagery of medium and low spatial resolution. The mismatch points can be efficiently reduced by the filter method. At last, only the glacier displacement signal $offse{t}_{ice}$, residual long-wavelength signal error $offse{t}_{res}$, and topographic effect $offse{t}_{top}$ are left.

3.2. DEM-Assisted Topographic Effect Removal

After all of the above procedures, the ice-flow related signal seems to be obtained. However, the topographic effect cannot be ignored with a large spatial baseline or rugged terrain [37]. Unlike the Polar Regions, topography in mountainous regions usually is complex, and rugged terrain will result in an obvious topographic-related component, especially in the slant-range direction with a large geometry baseline. Thus, an extra processing step of topographic effect compensation is needed.

A SAR sensor usually sends electromagnetic wave at a fixed rate determined by Pulse Repetition Frequency (PRF), and the echo is sampled also at a fixed rate corresponding to Frequency of Sample (FS). For single look SAR image pixel $\left(l,p\right)$, its range acquisition time $\tau$ and azimuth time $t$ can be written as [38]:

$$t=t_0+\frac{1}{PRF}l$$

(2)

$$\tau ={\tau }_0+\frac{1}{FS}p$$

(3)

where l and p are the row number and column number respectively in the two dimensional SAR image. $t_0$ is the acquisition time of the first line. ${\tau }_0$ is the delay time for the first sample in slant rage direction.

For the sake of simplicity, we assume that the radar flying direction (azimuth direction) is perpendicular to the radar look direction (slant range direction). According to the SAR geometry equations, the offset of same point in the master and slave can be expressed as following equations [38].

$$\Delta l=\frac{PR{F}_s}{v_s}{\widehat{v}}_s•\mathrm{{\rm P}}-\frac{PR{F}_m}{v_m}{\widehat{v}}_m•(\mathrm{{\rm P}}-{\text{S}}_0)+bia{s}_{az}$$

(4)

$$\Delta p=\frac{2F{S}_m}{c}\left|\mathrm{{\rm P}}-{\widehat{v}}_s({\widehat{v}}_s•\mathrm{{\rm P}})\right|-\frac{2F{S}_s}{c}\left|\mathrm{{\rm P}}-{\text{S}}_0-{\widehat{v}}_m({\widehat{v}}_m•\mathrm{{\rm P}}-{\widehat{v}}_m•{\text{S}}_0)\right|-bia{s}_{rg}$$

(5)

where $R_m^{\left(az\right)}=PR{F}_m/v_m$ and $R_s^{\left(az\right)}=PR{F}_s/v_s$ represent the pixel spacing of azimuth direction in master and slave full-resolution SAR intensity image. c is a constant. $bia{s}_{az}$ and $bia{s}_{rg}$ are constants for the initial offset in azimuth and range direction. ${\widehat{v}}_m$ and ${\widehat{v}}_s$ are the velocity unit of sensor in master and slave orbit. $S_0$ stands for the offset between the master and slave orbits. $R_m^{\left(rg\right)}=c/F{S}_m$ and $R_s^{\left(rg\right)}=c/F{S}_s$ respectively represent the range pixel spacing of the master and slave full resolution SAR image.

Differential of Equation (4) gives the azimuth local additional shift corresponding to an elevation on the ground-point position [38].

$$\text{δ}\left(\Delta l\right)=\left(\frac{R_m^{(az)}}{R_s^{(az)}}{\widehat{v}}_s-{\widehat{v}}_m\right)•\frac{\text{δ}\mathrm{{\rm P}}}{R_m^{(az)}}$$

(6)

where $\Delta h=\delta \text{P}\sin \theta$ is the topographic (vertical) variation above or below the reference altitude and $\theta$ is the look angle. We further assume that flight tracks are coplanar and present a small horizontal crossing angle $\alpha$. ${\widehat{v}}_s-{\widehat{v}}_m$ is in the across-tracking plane and $\left|{\widehat{v}}_s-{\widehat{v}}_m\right|\approx \sin \alpha$. The $\delta \text{P}$ is selected in the across-track plane and perpendicular to the range direction, in order to analyze the effect of a topography error at constant range and azimuth positions. With this hypothesis, the above equation becomes [38]

$$\text{δ}\left(\Delta l\right)=\sin \alpha \cot \theta \frac{\Delta h}{R_m^{az}}$$

(7)

Differential of Equation (5) gives the range local additional shift corresponding to an elevation on the corresponding ground-point position [38]

$$\text{δ}\left(\Delta p\right)=\left(\frac{R_m^{\left(rg\right)}}{R_s^{\left(rg\right)}}{\widehat{\rho }}_s-{\widehat{\rho }}_m\right)•\frac{\text{δ}\mathrm{{\rm P}}}{R_m^{(rg)}}$$

(8)

where ${\widehat{\rho }}_s$ and ${\widehat{\rho }}_m$ are the unit vectors that identify the master and slave (slant) range directions, respectively. And where $\Delta h=\text{δP}\sin \theta$. When the range distance is large compared to the baseline, we can assume that ${\widehat{\rho }}_s-{\widehat{\rho }}_m$ is perpendicular to the range and that its amplitude is $\left|{\widehat{\rho }}_s-{\widehat{\rho }}_m\right|\approx \sin \Delta \theta \approx B_{\perp }∕{\rho }_s$. Thus, the corresponding value becomes [38]

$$\text{δ}\left(\Delta p\right)\approx -\frac{B_{\perp }}{\text{ρ}\sin \theta }\frac{\Delta h}{R_m^{(rg)}}$$

(9)

According to different mechanisms revealed by the Equations (7) and (9), we can accurately estimate and compensate the topographical effect in both range and azimuth directions. Furthermore, it can be approved that the topographical related offset in range is larger than that in azimuth direction, because the cross angle between master and slave orbit is usually small enough to diminish the topographical effect in azimuth direction [37,39]. Therefore, the topographic effect should be eliminated in the slant-range direction with the help of external auxiliary DEM data of the study area (Figure 2).

Figure 2. Glacier surface motion field of the study area. (a) Displacement field in the slant-range direction with topographic effect. (b) Contribution of the complex terrain of the study area in the slant-range direction. (c) Last displacement field in the slant-range direction without topographic effect or residual orbit error. (d) Displacement field in the azimuth direction without residual orbital error.

Figure 2. Glacier surface motion field of the study area. (a) Displacement field in the slant-range direction with topographic effect. (b) Contribution of the complex terrain of the study area in the slant-range direction. (c) Last displacement field in the slant-range direction without topographic effect or residual orbit error. (d) Displacement field in the azimuth direction without residual orbital error.

After the topographic effect removal, a small long-wavelength residual signal sometimes appeared in the offset displacement field, mainly caused by the inaccurate modeling parameters for satellite orbit and sensor attitude. Such residual errors could be easily eliminated with linear or polynomial regression. Then, a smooth filter was employed to eliminate the noise in the displacement result, based on the combined constraints on modulus threshold and flow direction [3].

3.3. Error Analysis

The error in the intensity-based pixel-tracking displacement field is mainly generated in operations of image coregistration, parameter estimation for polynomial regression and topographic compensation. The histogram corresponding to residual error is shown in Figure 3. In order to evaluate the accuracy of glacier velocity, the mean and standard deviation (STD) value of the residual displacement, respectively 0.462 m and 0.404 m in 46 days, were computed based on the reliable pixels in the ice-free region around the mountain glaciers. If without topographic compensation operation, the mean and STD value should be 1.612 m and 1.082 m in 46 days. Furthermore, the statistical characteristics of residual error in the ice-free region, 0.462 m in 46 days, can be considered as the accuracy of final results. Therefore, the mean movement on non-glacial areas is well below the amplitude correlation accuracy.

Figure 3. Histogram of residual offset values in the non-glacial area.

Figure 3. Histogram of residual offset values in the non-glacial area.

The measured velocity 0.68 m∙d⁻¹ at the lower stake and 0.7 m∙d⁻¹ at the upper stake, which was carried out with DGPS during August in 2009 by Astrid, et al. [33]. In the corresponding positions on glacier surface, the ice velocity obtained with ALOS/PALSAR in winter is 0.78 m∙d⁻¹ and 0.50 m∙d⁻¹ respectively. In the lower part, the velocity from ALOS/PALSAR in winter is faster than that observed by DGPS and TerraSAR-X imagery in August. It is possibly caused by the much more precipitation and increasing thrust from the ice mass in the upper stream. On the upper stake, the measured ice velocity in this study is smaller than that of obtained by the DGPS. It is mainly due to the cold temperature and the basal slipping in the upper part is weaker than that in August. After all, due to the different seasons of velocity data and the variable ice motion in different seasons of year, the difference between these results can be considered as the reasonable diversity. The average velocity along the profile near the Fedchenko Glacier Station is about 0.85 m∙d⁻¹ in our results, which is well matched with the results from the long-term observation along the cross profile near the Station and the velocity measurements from the TerraSAR-X imagery [33].

4. Results and Analysis

The surface velocity field was successfully estimated with the available SAR image pair covering the glacial area of interest in the eastern Pamirs of Central Asia. The final displacement result confirmed certain basic general patterns in the motion of the glaciers and provided a better understanding of the glacier motion. For better understanding of the glacier activity, an extracted glacier surface velocity map was superimposed on the terrain relief map in Figure 4. Most of the large glaciers presented in Figure 5 are flowing approximately northward.

Figure 4. Geocoded glacier surface flow velocity field superimposed on a shaded-relief map of the DEM. The number of 3650, 4300, 4600 and 5000 indicates the elevation of corresponding locations.

Figure 4. Geocoded glacier surface flow velocity field superimposed on a shaded-relief map of the DEM. The number of 3650, 4300, 4600 and 5000 indicates the elevation of corresponding locations.

The ice flow field almost covering Fedchenko Glacier was accurately obtained to reveal its motion characteristics, except for the terminus because of the limited coverage of SAR data. But the motion on all the other parts of this glacier were clearly reflected in the final glacier flow map. The overall velocity pattern on the main trunk of Fedchenko glacier in winter 2007 indicates two regions from 3650 to 4300 m elevation and from 4600 to 5000 m elevation, respectively with mean value of 0.63 m∙d⁻¹ and 0.47 m∙d⁻¹ which is higher than the rest part of glacier. It is almost the same with the result presented by Astrid [33]. In addition, there are many small glaciers that can be observed in the glacier flow map and their flow patterns are also clearly exhibited. Generally, for small glaciers, the motion becomes greater with increasing altitude. On the tongue of these small glaciers, for example No. 1–6, ice flow velocity usually is less than 0.09 m∙d⁻¹ and sometimes stagnant. For some large glaciers, such as Fedchenko and Grumm-grzhimaylo glacier, the relationship between ice motion and altitude variation is much more complex, which is obviously shown by the profiles in Figure 6. It can be seen that the large glacier flow pattern is much more distinct than that of the small glaciers. The velocity in the original accumulation zone of the large glacier is less than that in the middle and downstream parts. Figure 4 reveals the large difference between ice velocities on different sections of the glacier. It shows relatively large motion with large fluctuations in the ablation zone of Fedchenko Glacier, whereas some other small glaciers are less active and their velocities decrease with decreasing elevation on the same section. This is mainly caused by the topography in the tongue part and the volume of ice from the accumulation zone.

Figure 5. Ice flow direction superimposed on the colored glacier surface velocity. (a) Fedchenko Glacier (b) Nalivkina Glacier and No. 5 Glacier (c) Grumm-grzhimaylo Glacier.

Figure 5. Ice flow direction superimposed on the colored glacier surface velocity. (a) Fedchenko Glacier (b) Nalivkina Glacier and No. 5 Glacier (c) Grumm-grzhimaylo Glacier.

Because of the constraints from complex terrain and valley geometry in the study area, there are many small tributaries contributing their ice mass to the main trunk of the large glacier, which can be intuitively demonstrated by the final velocity map in Figure 5. The flow directions of ice are also seriously limited by the complex terrain on both sides, which is obviously reflected in both Figure 4 and Figure 5. With the help of the glacier flow distribution, we can directly identify the active tributaries and further analyze the relationship between them and the velocity variation of the glacier trunk. The changes in velocity cross the flow line reveals that ice near the side-wall of the surrounding mountain becomes slow or even motionless due to lateral and bottom friction. Generally, the maximum ice flow velocity across the glacier can be achieved around the central flowline from the side-walls [3]. Similar phenomena are also observed in other large, wide glaciers [24,25].

Fedchenko Glacier and Grumm-Grzhimaylo Glacier are the two largest glaciers in the study area with comparatively fast flow velocities along the valleys. The others are short in length and show relatively simple variation in surface motion. Grumm-Grzhimaylo glacier not only has many tiny tributaries in the accumulation zone, but also has two relatively long, narrow tributaries in its middle section. The ice mass flows into the glacier trunk and slightly increases glacier surface velocity again from the junction between two large tributaries and the trunk of Grumm-Grzhimaylo Glacier. But the steep terrain is the key important factor causing relatively large velocity variation on the glacier surface.

In the optical image, the No. 5 and Nalivkina glaciers can hardly be separated from each other because of the saturation phenomenon caused by the serious continuous snow cover (Figure 1). But the glacier surface motion map in Figure 5b obtained in this study gives us an optimal way to distinguish one from the other. Based on the velocity map, the ridge between these two glaciers, indicated by the solid black line in Figure 5b, can be approximately obtained along the middle line in the stable area or zero motion gradient area. In Figure 5b, the ice is moving in the different directions on both sides of this line. Thus, the glacier surface motion map could present a novel way to outline the different ice mass accumulation region. Sometimes it could also play an important role in the glacier development research, because it is useful in outlining and analyzing the ice mass contribution area in the ice cap, where the terrain is relatively flat and the ice accumulation area is hard to identify.

Trunk No. 1–5 are small sub-glaciers that contribute ice mass to the trunk of the large glaciers, especially Fedchenko Glacier. Furthermore, they also have a positive or negative effect on the velocity variations of Fedchenko Glacier to some extent by ice mass contribution. From the velocity field of Fedchenko Glacier, it is obviously deduced that the ice motion was different at various parts and did not show a direct linear relation with the increase in altitude (Figure 6a). The velocity pattern reveals that the glacier generally flowed faster in the ablation zone than in the accumulation zone. Furthermore, in the ablation zone, the glacier shows a strong undulating movement. The Fedchenko Glacier reaches its maximum motion value at about 25 km up from the glacier terminus. The velocity distribution at two large corners exhibits complex spatial characteristics corresponding to the ice dynamics. The maximum velocity is slightly out of the central flowline because of the combined action of inertia and the effect of the surrounding mountains, which make sure the glacier flows along the valley. Five valuable large tributaries contributed ice mass to the main trunk of Fedchenko Glacier. Nalivkina Glacier, the biggest tributary, transferred ice mass to the accumulation zone of Fedchenko Glacier. The Bivachny glacier no longer flows into Fedchenko Glacier because it is stagnant on its downstream part, which can be observed and confirmed in the last displacement field result in Figure 4. However, Figure 4 also shows that Bivachny’s tributaries are still active and their ice becomes stable once it reaches the trunk of Bivachny Glacier.

The intersection angle between tributary and glacier trunk is also highly affecting the trunk surface velocity at the junctions, which is indicated by capital letter A–F in both Figure 1 and Figure 6a. In Figure 5a, when the angle of ice from some tributaries flowing into the trunk of Fedchenko Glacier is large enough, the ice movement would be slowed slightly at the junction of tributary and trunk because of the mass contribution with partly adverse direction, which counteracted the ice motion of the trunk, such as at the junction between tributary No. 2 and the trunk. Otherwise, the velocity increased when the ice flow direction of the tributary is generally the same as the trunk of Fedchenko Glacier at the junction, like the intersection of tributary No. 3 and the trunk. That is to say, the direction of the tributary flowing into the trunk plays an important role in accelerating or decelerating ice flow velocity of Fedchenko Glacier. These phenomena also partly reveal the dynamic characteristics of Fedchenko Glacier and its tributaries.

In order to investigate and give a clear, quantitative analysis of the ice velocity, some profiles of velocity were obtained along the approximate central flowlines of the glaciers, whose locations are indicated by lines in Figure 1. For analyzing the relationship between velocity and terrain, the glacier surface elevations were also included in the profile map in Figure 6.

Figure 6. Profiles of glacier surface motion in 46-day and their corresponding elevation variations. (a) Fedchenko Glacier; (b) Grumm-Grzhimaylo Glacier; (c) Nalivkina Glacier; (d) Bivachny Glacier.

Figure 6. Profiles of glacier surface motion in 46-day and their corresponding elevation variations. (a) Fedchenko Glacier; (b) Grumm-Grzhimaylo Glacier; (c) Nalivkina Glacier; (d) Bivachny Glacier.

The velocity values and distance from the terminus are shown in Figure 6 along the central flowlines, starting from the ablation zone (glacier terminus) up to the accumulation zone of the glaciers. Although the glacier surface elevation continually decreased down-valley, the derived velocity profiles display significant variations along the glaciers’ central flowlines, especially in the ablation zone (in the lower sections) (Figure 6a). Along the central flowlines of the glaciers, there were multiple local maxima and minima, indicating that the motion of ice mass in relatively large glaciers with long length is much more complex than that of the relatively small and short glacier. Generally, the velocity distribution pattern of large scale glacier is much more variable than that of the small scale glacier in mountain areas, which is mainly caused by the rugged terrain around the glacier and various supplement of ice mass. The mean velocities of most glaciers in study area fell between 0.08 m∙d⁻¹ and 0.57 m∙d⁻¹, with ice showing a maximum velocity of 0.85 m∙d⁻¹ in Fedchenko Glacier.

Both Figure 4 and Figure 6 indicate that the longer the glacier extends in mountain area, the more complex the corresponding surface velocity is, which mainly due to the rugged terrain experienced by the moving ice. This is obviously demonstrated by Fedchenko Glacier in mountain area, the longest glacier outside of the polar region. Both Figure 4 and Figure 5a show that the velocity distribution was much more complex on Fedchenko Glacier than that of the smaller glaciers. Besides the reasons mentioned above, a few large tributaries imposed their positive or negative influence on the velocity of Fedchenko Glacier. Moreover, both the velocity field map and the longitudinal profiles show that the velocity in the upper part was less than that in the middle and downstream section of the glacier, mainly due to the relatively flat topography and colder temperature in the accumulation area.

Additionally, the complex and rugged terrain in the study area not only makes the glacier shapes and their surface velocity distributions vary significantly, but they also introduces a serious topographic effect in the velocity field extracted based on the SAR data, which has a negative influence on the accuracy of the derived ice velocity field. According to the basic principal hypothesis of the pixel-tracking method, the topographic effect should be rightly and precisely compensated for. Otherwise the accuracy of glacier surface velocity field will be seriously affected. In this paper, an external SRTM DEM version 4 with 90-m resolution was used to efficiently remove the topographic-related displacement, and then the residual offset in the free glacial region could be used in the accuracy estimation.

5. Discussion

The equilibrium-line altitude (ELA) of a glacier is sensitive to climate change and tightly correlated with solid precipitation and air temperature. According to the theory of glacier kinematics, glacial surface velocity usually reach the maximum near the ELA [40]. Therefore, the elevation where maximum glacial surface velocity occurs will, to some extent, indicates the location of the ELA. The accumulation area and the ablation area, revealed by the glacier motion distribution, give us an opportunity to understand the situation of the whole glacier. The long glacier, such as Fedchenko glacier, should have the large accumulation area or accumulation–area ratio (AAR), which is necessary in keeping the glacier health.

Unfortunately, the rugged terrain around the mountain glacier usually presents the topographical related value in the results of pixel-tracking method. In order to extract an accurate ice velocity, the effect of elevation variation in the pixel-tracking method was analyzed and estimated with the help of external DEM data in study area. Due to the small cross angle between master and slave SAR orbits, the rugged terrain usually presents the weaker effect in azimuth direction than that in the range direction. Therefore, the topographical effect compensation should be always done in the range direction in complex topographical area, not to the azimuth direction. Furthermore, the topographical effect in range direction also can be ignored when the spatial perpendicular baseline is small enough.

Along the glacier flow direction, the ice velocity is changing heavily. It is because that both the rugged terrain and the ice mass from the tributaries in the accumulation part make the motion of ice diversity in large scale glacier. Therefore, both variable gradient associated with the rugged terrain and gravity of ice mass, which are the important influence factors for ice motion, make the glacier surface flow distribution much more complex. Generally, the low gradient with flat terrain would make glacier flow slowly, especially in the origin of the glacier. The friction from the surrounding mountain also influence the ice flow variation across the glacier, which usually make the ice reach its max velocity value near center part along the cross profile.

The variation of velocity on glacier surface generally illustrates the dynamic of glacier, which could be used for predicting the ice motion related hazard. Although there is an observatory station for glacier studies around the middle section of Fedchenko glacier, little is known about the overall behavior of the glaciers in this area because of their huge coverage and long length, which make field investigations performed in this glaciated area suffer heavily from the harsh weather and complex terrain conditions. Only studies based on the remote sensing method, rather than ground measurement, can provide ice movement information covering the whole area of interest.

In the eastern Pamirs, the glaciers are almost covered by snow, especially in the cold season because of a large amount of precipitation. In addition, the study area has a high average altitude and is often covered by clouds. All of this makes optical imagery unsuitable for distinguishing the contours of glaciers, not to mention the motion pattern of ice. ALOS/PALSAR imagery, acquired with an L-band operational SAR sensor, is hardly influenced by the cloud and snow cover and became the only useful data source for investigating the glacier dynamics. Compared with C- and X-band, L-band can provide much higher quality data in snow-covered glacial area because of its deep penetration. The vertical ice motion profile indicates that only velocity near the bed, rather than near the glacier surface, experiences the most significant changes [35]. Therefore, the differences in the penetration depths corresponding to the different wavelengths will not affect the derived glacier surface motions except for the ability in maintaining the high enough temporal coherence. In mountainous glacial areas, the complex, rugged topography presents trouble for accurately extracting movement on glacier surfaces. But it can be efficiently modeled and removed from a glacier’s flow field with the assistance of external SRTM DEM for accurate estimation of ice displacement field.

Glacier retreat is pronounced in the study area [31]. Glaciers are expected to be especially sensitive to present-day global warming and climate change, and as glacier size decreases, glacier velocity also gradually decreases. This phenomenon is significant in study glaciated region. This trend is obviously not only in the eastern Pamirs but also in most glaciated areas in the Himalayas. Owing to glaciers’ remote geographic location and harsh climate conditions, SAR images will be an increasingly important glaciological tool in both the extraction of ice surface velocities and recognition and explanation of ice motion events. Pixel-tracking is a useful method for obtaining glacier velocity measurements using SAR imagery, especially under difficult geographic conditions [4,22]. The integration of SAR remote-sensing and regular in situ measurements of representative glaciers in the eastern Pamirs is therefore highly valuable and recommended for future glaciological research in the context of climate warming [3].

The high accuracy of glacier surface velocity also plays a key role in the study of glacier dynamics and provides a potential opportunity to understand the activity of the ice, especially in remote regions. In addition, it is apparent that the glacier boundaries can be clearly identified in the glacier velocity map. It is helpful to outline the glaciers always covered by the snow, which cannot be easily distinguished from the surrounding snow-covered mountains, especially in optical remote sensing data.

Furthermore, the velocity distribution map to some extent can also assist us in locating the regions with a high possibility for crevasses to occur. Crevasses mainly form in areas with tensile stress, usually corresponding to a part where velocity is greatly increasing. Besides, it provides useful information for field investigations because the identification of potential dangers on the glacier surface can help us optimize the routine to avoid concentrations of crevasses in expeditions to the study area. In addition, there is much more work should be carried out for comprehensively understanding the long-term variation of the glacier movement in study area.

6. Conclusions

This paper presents a method to efficiently provide an accurate ice velocity distribution map for several glaciers located in the eastern Pamirs of Central Asia using pixel-tracking methods assisted with external SRTM DEM version 4 data for compensating for topographical effects, based on ALOS/PALSAR images acquired on 2 January and 17 February 2007. Given the relatively simple and efficient procedures and the excellent accuracy of the derived results, the pixel-tracking method with SAR imagery identified in this study is demonstrated to be a suitable approach for estimating glacier surface velocities, especially in cloud- and snow-covered glaciated mountainous areas where visible optical imagery is problematic. The methodology presented in this study, including DEM-assisted topographic effect removal, can be readily applied to accurately extract glacier flow velocity distribution in other mountainous glaciated areas, given the availability of SAR remote sensing data. Thus, the method can be employed to accurately monitor mountainous glacier velocity routinely in remote, hard-to-reach areas with rugged terrain.

As shown in the obtained velocity field, both the spatial distribution pattern and velocity variations of the ice flow reveal that glacier velocity is seriously confined by surrounding valley geometry and local topography in mountainous terrain. The complex, rugged terrain not only limits the glacier flow patterns, but also presents a considerable topographic effect that can be efficiently compensated for with the help of both external DEM data and the geometric relation between two SAR acquisitions. Besides the terrain and climatic factors, surface velocity variation is also lightly affected by the ice mass flow from tributaries. They have a positive or negative effect on the ice velocity depending on the angle between trunk and tributaries, which makes the ice velocity more complex in the intersections of trunk and tributaries than other sections (see Figure 4 and Figure 5). In order to detect seasonal fluctuations of the glacier surface motion and to further estimate the mass balance, a large amount of SAR imagery covering the same area of interest is required to be routinely obtained for comprehensive investigation of the surface velocity. Furthermore, regular, long-term monitoring of alpine glacier surface velocity in Central Asia is vitally necessary to provide valuable sensitive indicators of current and future climate change for environmental analysis.