Glacier surface velocity estimation in the West Kunlun Mountain range from L-band ALOS/PALSAR images using modified synthetic aperture radar offset-tracking procedure

Abstract Glacier movement is closely related to changes in climatic, hydrological, and geological factors. However, detecting glacier surface flow velocity with conventional ground surveys is challenging. Remote sensing techniques, especially synthetic aperture radar (SAR), provide regular observations covering larger-scale glacier regions. Glacier surface flow velocity in the West Kunlun Mountains using modified offset-tracking techniques based on ALOS/PALSAR images is estimated. Three maps of glacier flow velocity for the period 2007 to 2010 are derived from procedures of offset detection using cross correlation in the Fourier domain and global offset elimination of thin plate smooth splines. Our results indicate that, on average, winter glacier motion on the North Slope is 1     cm / day faster than on the South Slope—a result which corresponds well with the local topography. The performance of our method as regards the reliability of extracted displacements and the robustness of this algorithm are discussed. The SAR-based offset tracking is proven to be reliable and robust, making it possible to investigate comprehensive glacier movement and its response mechanism to environmental change.


Introduction
Glacier change is an important indicator of global climate change, especially for mountain glaciers, which are highly sensitive to variations in temperature and precipitation. [1][2][3][4] Because most glaciers are located in sparsely populated locations far away from major human activities, the relationship and response mechanism between glacier movement and environmental change can be more relied upon in the study of global change issues without human influence. In addition, the climate change occurring at present is shifting the geomorphodynamic equilibrium of glaciers and intensifying glacier-related natural disasters, such as landslides, flood, and permafrost creep. 5,6 Therefore, it is crucial to detect and monitor glacier movement in many areas. Most glaciers are located in high-altitude mountain areas, and it is difficult to measure large-scale glacial variables through conventional ground surveys. Up to now, only 1% of the world's existing temperate glaciers have been monitored, mostly by ground measurements, which often provide information only once or twice a year at a few points. 7 Remote sensing techniques have the potential to obtain regular observations of glacier movements. In particular, synthetic aperture radar (SAR) data have been efficiently applied to detect glacier movement regardless of weather conditions. 8,9 Current motion detection based on the SAR images contains two well-tested approaches that use different information contained in SAR data-phase and intensity. One is differential SAR interferometry (D-InSAR), in which two complex SAR images acquired from slightly different orbit configurations and at different times can be combined to exploit the phase difference between the signals. 10,11 Because the interferometric phase is related to the topography and surface displacement of the image pairs, the mean glacier velocity along the slant-range direction can be extracted after removing the topographic phase. Several studies have demonstrated that differential interferometry is suitable for obtaining the whole-glacier velocity field quickly. 12,13 However, although the D-InSAR procedure chain is built from standard InSAR algorithms, it often fails to acquire useful fringes due to the decorrelation caused by complex glacier topography, phase noise, atmospheric heterogeneities, or large glacier deformations. 14 Another technique is SAR image offset tracking, also known as the cross-correlation optimization procedure, which calculates the displacement of identical objects from two scenes with different imaging times. The SAR offset tracking is based on the cross-correlation algorithm, which is a similarity measurement used in image registration for determining homologous points and detecting the offset. 15 There have been a few studies which apply the offset-tracking method to SAR or optical satellite imagery for large-scale monitoring of mountain glaciers, and their results validate the efficacy of this novel approach. [16][17][18][19][20] Compared to the precision of the deformation derived from D-InSAR (millimeter level), the offset-tracking method provides less precise measurements on the order of centimeters. 21 However, the accuracy is still within 1∕5 of the size of a pixel, which is acceptable for the estimate of rapid flow areas where there are no available ground survey data or interferogram fringes that could feasibly be interpreted. 22 In addition, the offset-tracking procedure can provide displacements in the two-dimensions of the range and azimuth directions. The D-InSAR method yields the offset in the slant-range direction only; therefore, the SAR-based offset tacking is more suitable for producing comprehensive twodimensional velocity fields. 23 The correlations between two images contain not only offsets related to the movement of ground objects but also several factors independent of geographical phenomena. The attitude of imaging sensors and local topography is considered to yield major distortions, and these offsets generally show wave artifacts in the across-track direction of the image. The step where residual offsets are removed involves modeling the unwanted components using attitude parameters and digital elevation model (DEM) information, or using a polynomial fitting function. However, the geometric distortion of the SAR images in mountain areas is quite complex, due to uncertain imaging parameters or topographic relief effects, and, therefore, it is difficult to ensure the quality of removing the additional offsets from real glacier motions using these approaches for simulating global offsets.
In our study site-the West Kunlun Mountain range-glaciers are located above the plateau surface of 4500 m and generally reach altitudes of more than 6000 m. 24 High-altitude topography reduces the accessibility of routine ground measurements and the reasonable interpretation of valley glaciers using interferograms. In addition, this glacier area lacks a glaciological investigation base using ground measurements or satellite data. Except for the large-scale expedition to the West Kunlun mountain glacier area, which aimed to study glacial parameters and changing behavior in the 1990s, there have been no other relevant intensive studies up to now. In this article, we use the SAR intensity offset-tracking method to derive the glacier surface velocity field in this region. Meanwhile, we provide a modified offset-correction method by introducing another fitting algorithm and conduct an evaluation of the method's precision and applicability for monitoring glacier dynamics.

Study Area and Datasets
The West Kunlun Mountain range is one of the densest areas of glaciers in China, and there is a dense distribution of glaciers along the mountain ridge. The total glacier area around the Kunlun Peak (7167 m) is about 5640 km 2 . 25,26 The West Kunlun Mountain range is located on the border between Xinjiang and Tibet in the west of China (34.5°to 36°N and 79.8°to 82°E). According to the report "Global Land Ice Measurements from Space (GLIMS)," there are 3165 glaciers in this area; among these, the Duofeng glacier is the largest with a length of 26.8 km and an area of 251.7 km 2 . 27 This study focuses on the eastern portion of the West Kunlun Mountain range, which includes six well-known glaciers and two ice caps [see Fig. 1(a)].
The ALOS PALSAR imagery used in this study is illustrated in Fig. 1(b). In this imagery, the overall shape and detailed characteristics of the glaciers, such as debris and glacier tongues, are well demonstrated, which makes the data suitable for velocity extraction using homogeneity matching in the offset-tracking method. The eight marked glaciers and ice caps were chosen for velocity detection based on the ALOS PALSAR data acquired in winter during the period 2007 to 2010 (Table 1) to investigate local glacier changes and to explore the tools available for glacier observation.

Modified SAR Offset-Tracking Procedures
We improved the conventional SAR imagery offset-tracking method by choosing the search window of 64 × 64 pixels suitable for our study site and by fitting a global offset using the thin-plate spline (TPS) algorithm. We obtained the offsets for pixels that had moved between the acquisition of the images in the ALOS PALSAR pairs. The glacier surface velocities could be extracted from the measured offsets and the temporal separation. The glacier flows in winter for the period 2007 to 2010 were derived for the eastern section of the West Kunlun Mountains.   Figure 2 illustrates the flow chart for acquiring glacier surface velocities in two horizontal directions. After the global offset estimation for the image pairs, the two images could be coarsely aligned in the same coordinate system. The second step was to choose a proper size of the patch containing sufficient features in one image (called the master image) and then search for its corresponding patch in the slave image, in order to obtain a subpixel image motion vector using the upsampled cross-correlation method in the Fourier domain. 28 The crosscorrelation algorithm is well established for matching similar features but the offset originally produced contains global misregistration with a general trend in both the range and azimuth directions. 29 The intermediate processed images shown in the flowchart (Fig. 2) indicate that the original offsets in the two directions have an obvious striped component distributed over the whole image. After modeling the global offset with a fitting algorithm suitable for the study of mountain glaciers and appropriate for the SAR imagery, relatively accurate displacements can be derived. (These images are shown in the last step of the flow chart.)

Modified Global Offset-Correction Method
The horizontal offsets acquired in step 3 contain additive distortions caused by the imaging sensor's attitude and the topographic complexity and usually require two corrections: an orbital correction based on fitting a bivariate polynomial function to offset fields for an area of stable ground and a topographic correction based on a simulation produced by a DEM. Since it is difficult to quantify the distribution of each factor, the corrections are often simplified to one step by simulating global offsets using a bivariate polynomial function or a binary quadratic polynomial function, under the assumption that no displacement occurred in most parts of the image. However, this process is often time consuming. In addition, these functions for modeling rigid deformation can explain complicated surface features relatively simply.
TPS refers to the physical analogy involving the bending of a thin sheet of metal, and this is becoming an efficient algorithm for curved surface interpolation based on discrete points. TPS can be used to build both rigid and nonrigid transformations. It has been widely used for complex surface modeling in image alignment, surveying and mapping, or for climatological and hydrological interpolation. [30][31][32] Compared to the polynomial functions, TPS can express complex deformations more smoothly and elastically. In addition, TPS can be applied over whole images containing moving regions, thus reducing processing time.  Figure 3 illustrates the initial offsets and the global offset simulated by TPS, and their trends in both the range and azimuth directions. The simulated global offsets have the same distribution trend as the initial offsets without accurate registration, but also demonstrate the local fluctuations. Specifically, the two offsets have the same smooth ascending trend within [−2 AE 0.5, 2 AE 0.5] in the range direction and have a similar trend line-not changing by more than 0.25 pixels in profile-along the azimuth direction. In addition, TPS is able to obtain specific features, for example, in the region of azimuth 100 to 200 and range 1 to 20, and the two offsets exhibit the same fluctuation if the oversized values, which are regarded as noise, are ignored. Therefore, in our study, this modified procedure is used to remove both systematic errors and the influences of complex topography.

Accuracy Assessment
The accuracy of the offset-tracking method for the SAR imagery is limited by several factors, such as decorrelation due to long time intervals, the small matching patch size, or objects being unsuitable for registration. 33 Errors caused by these factors could be eliminated through detailed analysis and by finding the acceptable threshold values of the time separation or patch size by trial and error. Also, objects such as water bodies in an image could be rejected because of the inconsistent backscattered signal at the different acquisition times. However, some factors could not be eliminated due to their random nature and complexity-these include the noise signal and distortion in the imagery. It is impossible to provide a theoretical bound for this error but alternatively, the accuracy assessment in this study aims to estimate glacier offset values by comparing the displacement values generated over stationary areas. 34 There is a lack of large-scale accurate ground measurements at the test site in the past 10 years and thus the results will be evaluated empirically by considering the potential motion according to geological conditions. . The results are given as the Euclidean distance of the offsets (in pixels) measured in the slant-range and azimuth directions. In this area, glaciers are distributed along the north-south direction, and the glacier motion field is generally along their orientation, in accordance with glacial geology. Although the offsets in the slant-range direction are not >1 pixel on average and <2 pixels in fast-flowing areas, this component is a minor one in the glacier velocity field compared to the main flow-line component in the azimuth direction. In addition, based on ALOS PALSAR imaging theory, the range resolution on the ground is about 7.5 to 8.3 m∕pixel compared to 3.5 m∕pixel in azimuth, so that the range component carries more weight in the calculation of the velocity. Therefore, in order to reduce the error caused by the quality of the SAR images, we eliminated the offset components along the azimuth, which were several times larger than the ones along the range direction.

Velocity Map of Glacier Flow
We have listed detailed information for six glaciers and two ice caps measured by the offsettracking method in Table 2. We found that: 1. Winter displacement of glaciers on the North Slope is larger and more consistent than that on the South Slope and the velocities of glaciers at similar longitudes show a reduction of about 1 to 2 cm∕day from north to south. The reason for this discrepancy is probably that the glaciers on the northern slope mostly have longer trunks corresponding to the greater difference in elevation than the ones in the south. For example, the Northern Slope Glacier has an elevation that ranges from 6300 m in the glacial accumulation part to 4500 m at its front. This can be compared with the Chongce Glacier, which has a southerly aspect and an elevation ranging from about 6300 to 5300 m. To facilitate comparison of the opposite aspects, we extract the elevation profiles of typical glaciers on different slopes (Fig. 5). The velocities of glaciers are related to the changes in elevation according to the linear regression analysis: the elevation declines more on the north with an average slope of the regression of −0.93, and for the glaciers with slow movements the topographic changes are gentler with a slope of −0.43. 2. Duofeng glacier on the North Slope has the highest forward velocity, which is about 10 cm∕day on average. The lowest velocity is for the Yulong glacier in the northeast, about 1.2 cm∕day, and there was no recognizable flow field extracted along its trunk. Moreover, Chongce Glacier has an obvious flow field in the upper part of the main trunk only, and the velocity of the glacier tongue is very close to that of stationary land. These results are partly related to the elevation gradient of the glaciers-steeper topography means a greater gravitational field acts on the glacier body and produces the glacier's sliding motion. These three glaciers-Chongce, Yulong, and Duofeng-have similar lengths of 28. 7, 30.9, and 31.0 km, respectively, but the Duofeng glacier has an elevation range of nearly 1000 m compared with the 400-to 500-m altitude difference of the two glaciers that do not have a fast flow field. 3. Two large ice caps in this region have the lowest velocity without showing any regional characteristics. The survey of the Guliya ice cap in 1991 indicated that it was composed of three arched flat ice caps and several ice tongues extending from west to east. 35 According to glacier dynamics, different types of ice will demonstrate their own movement features; this is especially true for ice tongues, where there will be relatively obvious flows along their main trunks. Although no ground survey of the Chongce ice cap has been made, from the Landsat thematic mapper (TM) images it is obvious that this ice cap also contains several glacier tongues. The extracted displacements perhaps contain processing errors due to the topographic conditions and the range in the amount of backscatter from flat snow and ice surfaces and require more research and observation for validation.

Duofeng glacier is a branched valley glacier located on the North Slope of the West Kunlun
Mountains. It is a large glacier with a length of 31 km and an area of 251.7 km 2 . 36 The surface velocity field of Duofeng glacier was extracted for its main trunk along the direction of its orientation (Fig. 6). In order to find the influence of changes in glacier elevation, we compared the displacements along a longitudinal profile of the three image pairs (Fig. 7). The main trunk of Duofeng glacier has a velocity field that averages 13 cm∕day in winter (December to January). The motion field appears to have a higher velocity in the upper elevations of the glacier body, with a continuous reduction in velocity to the front of glacier tongue. Figure 6 shows a similar trend of glacier flow during the period 2007 to 2010-the velocity is relatively stable at around 10 to 15 cm∕day along the main flow line with a decrease around an altitude of 5500 m and an increase around an altitude of 5300 to 5400 m. The higher velocity field in the upper section just below the snowline of 5760 m is forced by the gravity field of thick ice and snow in the accumulation areas above the snowline. Another high-speed region is potentially affected by the four branches sliding down from the adjacent sides of the mountain. Comparing

Guozha glacier
Guozha glacier is on the South Slope and has an elevation of 5390 to 6530 m, covering an area of 33.47 km 2 . 31 From the comparison of horizontal displacements (Fig. 8) and longitudinal profiles (Fig. 9), we found that Guozha glacier moves an average of 9 cm∕day in the period December to January. The velocities above 5800 m are much slower in all three profiles because these are the accumulation areas, where the major glacier movements are the accumulation of snow layers. Velocities increase to greater than 15 cm∕day around 5600 to 5800 m, probably as a result of the sharp elevation change and because of a branch flowing into Guozha from the Chongce ice cap to the east.

Result Assessment
Because of the restrictions due to the geographic conditions, insufficient comprehensive ground measurements or research by remote sensing observations have been conducted in the past  10 years. It is, therefore, difficult to carry out an accurate validation of the displacements obtained by the SAR offset-tracking methods. We evaluated our results using two criteria associated with the feasibility of the glacier velocity results and the techniques used.

Residuals in the stationary areas
The accuracy of the velocities can be first assessed by a null test over motionless, ice-free areas of the images. 22 For well-registered image pairs, the offsets over stationary areas are meant to be as small as possible. We masked the glaciated areas in the SAR images in order to conduct the residual analysis on the areas that were assumed to be stationary. In the SAR images, the expression of regions covered by snow and ice is too complicated to distinguish the glaciers' borders explicitly from the snow, especially in mountainous areas with rough topography. To avoid misclassifying the stationary areas, the mask used covered all the snow and glaciers in the mountain range as mapped by the Landsat TM data, which performs better in imaging natural objects. The residual statistics are listed in Table 3, including the mean values, the standard deviation and the ranges of the major offsets, which were collected from the histograms of the offsets in the three image pairs. We found that the offset residuals have a range of around 0.11 to 0.16 pixels, which means that the average error in the stationary areas could be under 1∕5 of a pixel over 46 days. The possible residuals from the offset-tracking method will not be more than 1∕4 of a pixel for a period of 46 days. Although there were some residuals left after removing the global offset, the value was within an acceptable range compared to the large extent of the glacier and the research area, taking account of the complexity of ALOS PALSAR imaging in mountainous places, which makes it hard to register images accurately.

Consistency of glacier velocity field extracted from the overlap
Assuming that the method applied to glacier areas in the SAR data is robust enough, an identical algorithm using different parameters must be able to produce relatively consistent results for the same areas. We divided the whole study site into two parts to be processed separately, with slightly different steps for the searching window [ Fig. 10(a)]. Because imaging attributes, geographical conditions, and objects in motion were the same in the overlap region, the results of the SAR-based offset-tracking method were affected mainly by the applicability of the algorithm to glacier observations based on the SAR data. We extracted one glacier flow field within the overlap to compute the overall offset and longitudinal profile offsets [Figs. 10(b) and 11]. The mean offsets differed by 0.25 pixels and the difference between offset values for the majority of the area was <0.5 pixels. The two longitudinal profiles show that this procedure could provide  Fig. 11(b). Fig. 11 Quantitative analysis of glacier flow offset in Fig. 10(b).
observations of glacier motion patterns using the trends that are the same in these two profiles; however, there is a deviation that reaches 0.5 pixels, which is likely caused by the constraints of ALOS PALSAR imagery and the limits of the interpretation of complex glaciers simply from intensity information in the SAR data. It means that the potential deviation using this method based on ALOS PALSAR data is <0.25 pixels (1.9 cm∕day) on average and that the maximum on the center line of the glacier flow is under 0.5 pixels (3.8 cm∕day).

Discussion
The study area is ideal for the processing and validation of SAR intensity offset tracking for glacier monitoring because of (1) the diversity and mass distribution of the glaciers, (2) the visibility and good presentation of glacier information observed from ALOS PALSAR, and (3) the high-relief topography, which limits the application of the SAR interferometry and dense ground measurements. Few published ground observations were made or detailed research of glaciers in the west Kunlun mountain range conducted in the last 10 years and, therefore, our results could complement glacier studies at this site and provide a reliable tool for such applications. Our results indicate that glaciers on the North Slope are moving faster than those on the South Slope. These results are supported by previous geological investigations which show that the North Slope of the Western Kunlun Mountains are severely cut and highly shadowed. 37 Because of the high altitude and high latitude, the temperature of the environment and of the glacier ice is very low and is not conducive to glacier motion. The glacier motion, therefore, is highly affected by the gravitational field produced by the different topographical conditions on the two slopes. 38 However, the specific velocity values obtained cannot be confirmed by previous measurements, which were either made 20 years ago or covered only a small region. Our study provides a reliable method of extracting regular and large-scale glacier velocity fields that are reasonable and also stable enough for further analysis to the same standards. In order to meet the requirements of the relative variable analysis, the performance of the SAR-based offset-tracking method must be evaluated in detail. We chose thin-plate smoothing splines to remove the global offset, taking advantage of its effectivness in modeling the gentle deformation by discrete points. This is different from conventional procedures, such as registering images using ground control points or imaging geometry, or removing a global offset using a polynomial expression. The offset maps of the three image pairs indicate that the global offset was effectively removed as displacements in glacier areas can be easily observed and there is no illogical regional offset distributed over the maps. By looking at the statistical distribution of the offset values, it can be seen that 80% to 90% of the offset measurements for stationary pixels are within the range 0.00 to 0.16 pixels and about 10% are <0.25 pixels, which means that the residual error will be not >1.90 cm∕day. The possible causes of the residual error may be random signal noise, a physical change in the entity observed during the time interval between the acquisition of the images or the complex fraction of the global offset, which is hard to extract. However, compared to the ALOS PALSAR imagery pixel size of 3.5 m in azimuth and 7.5 to 8.3 m in ground range, the accuracy of this procedure is acceptable for the requirements of further analysis.
Another evaluation indicator is the assessment of the robustness of this method, which is concerned with whether it is suited to large-scale observations over a long period. This uncertainty in the results is caused by image matching ambiguity, which is difficult to quantify. 39 Effective feature matching is affected by whether the features in the searching window can correctly reflect the essential characteristics for matching, or whether the characteristics could be interpreted by mathematical expressions. The process of coherence feature matching remains ambiguous due to the size and the sampling step of the search window and the ability of algorithms. In our research, we propose a complement to the evaluation of coherence matching based on the SAR features by reprocessing the same areas with different imaging times or with different searching window steps.
1. The data in our research, three image pairs for the period 2007 to 2010, provide the comparisons between different times. The three displacement results, both between winters and for times 1 year apart, were expected to barely differ from one another due to the short time differences and the slight variations in temperature and precipitation. Therefore, if these results give similar displacements, the algorithms can be considered as capable of calculating real and reliable values of the offset from the SAR matching coherence affected by different imaging parameters and signal noise. As shown in this article, the three results are close to each other over the whole offset map in terms of the statistics for the glaciers and stationary areas, and in the change in the local glacier flows. 2. Comparing the overlap of the same processed data with different search steps the deviation of the method can be found, meaning a quantified ambiguity of glacier feature matching applied to the ALOS PALSAR data. As discussed in Sec. 4.3, the results have similar characteristics with an average deviation of 0.24 pixels in 49 days (1.8 cm∕day), reflecting the potential error in the comparative study of glacier motion.
Although the results could not be validated numerically, the method described in this article is able to provide useful displacement values, which are logically consistent with geological and glaciological theories. Our results confirm that this method is suitable for use in glacier research, providing credible displacement distribution with acceptable errors. In order to optimize this method and adapt it to the characteristics of both glacier and the SAR data, more ground measurements for validation are needed in further studies.

Conclusions
We presented an initial study that used the SAR-based offset tracking for extracting glacier motion in the West Kunlun Mountain range from three pairs of the ALOS PALSAR images to derive the glacier flow velocity and to evaluate the capability of this method. This research indicates that higher flow velocities occurred on the North Slope of the mountain range than the south, confirming that topography is the main factor in glacier motion. The error in the results is considered to be <1.90 cm∕day and the error in the comparisons of glacier motion is considered to be 1.8 cm∕day on average. Therefore, our methodology can be applied in glacier research without ground control points and has the potential to provide glacier change analysis over a long period at a large scale. Shiyong Yan is a PhD student at the Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences. He received his BS and MS degrees in remote sensing technology respectively from the Chang'an University and Hebei University of Engineering in 2006 and 2009. His research interests include digital image processing, InSAR and PS-InSAR technology.