The south inylchek glacier activity analysis over a ten-year interval in surface velocity with multi-source spaceborne imagery

Glacier movement is one of the most important characteristics in describing mountain glacier activity, which is a sensitive natural indicator of climate change. However, the short-term glacier movement with a single data source could not precisely and sufficiently demonstrate the response of glaciers to climate change. In order to extract the reliable signal corresponding to climate change, the long-term monitoring of glacier movement should be widely exploited. This paper presents the glacier movement distribution of the South Inylchek Glacier by the improving pixel tracking algorithm with both synthetic aperture radar and optical satellite images in 2007–2008 and 2017–2018. The analysis in spatiotemporal characteristics of the glacier velocity indicates that the South Inylchek Glacier remained almost stable in a ten-year interval with the average velocity of 32 cm/d, which are computed with quasi-synchronization multi-source imagery. And the consistency of the velocity results was also verified with both optical and SAR imagery during the same period. Therefore, it would be valuable in expanding the research temporal cycle of glacier movement by uniting multi-source data. And the long-term glacier movement should be exploited for further analysis in mass balance prediction and assessing the climate change effects in the High Mountain Asia.


Introduction
With the increase of extreme weather conditions, as a natural indicator of global climate, glaciers in low-latitude and high-altitude regions have widely gotten people's attention. Glacier movement monitoring as well as corresponding variation could not only reflect material balance variation and surge, but also acts as a predictor for glacier disasters. Besides, as one of the important characteristics, glacier movement is closely related to its mass balance and stability. Furthermore, the long-term variation of glacier movement around the world would provide better insight into the spatiotemporal evolution of glaciers and global climate change (Paul et al. 2015;Shangguan et al. 2015;Tong et al. 2018).
Currently, most studies of glacier movement are usually revealed with a single data source (Paul et al. 2017;Sun et al. 2017;Shukla and Garg 2020). And there are few glacier movement studies based on the combination of multi-source remote sensing images. Generally, the optical imagery with better storage has been widely applied to extract glacier velocity in polar and high-altitude regions (Berthier et al. 2007;Nobakht et al. 2014). Despite its coverage scale and storage capacity, optical imagery is often hampered by inclement weather conditions and the mountain glacier surface is partly covered by snow all year round , which easily leads to an oversaturation phenomenon and destroys the integrity of results (Yan et al. 2019). Therefore, SAR imagery was employed in surface displacements estimation, thanks to more launches of radar satellites. And SAR imagery with high spatial resolution provides all-time, all-weather services, free from any cloud cover. Besides, SAR imagery with high spatial resolution could be exploited by differential interferometric SAR (D-InSAR), multi-aperture interferometry (MAI), and pixel tracking algorithm(PT) for more accurate glacier flow estimation 1 3 749 Page 2 of 12 with the phase and intensity information of SAR images and the external data (Strozzi et al. 2002;Yan et al. 2016;Tong et al. 2018). Spaceborne SAR then has been widely used for the multi-temporal analysis of glacier movements, yet such existing problems as geometric distortions, the oblique projection and the complex terrain condition would severely limit its application (Berthier et al. 2005;Shah et al. 2019).
Generally, the longer the observation temporal period, the more accurate the results revealing the characteristics of climate change. However, it is difficult to meet the requirements of such research with a single data source because of the limitations of satellite operation and poor image capability. There are few comparative analysis with a temporal interval of more than 10 years in previous studies (Huang and Li 2011;Paul et al. 2017). Nowadays, with the development of remote sensing techniques, it is possible to carry out glacier movement estimation with multi-source imagery for longterm glacier investigation (Lv et al. 2020;Peng et al. 2021;Singh et al. 2021). Thus, both the optical and SAR data with nearly a year interval from 2007-2008 to 2017-2018 are selected respectively and applied for glacier movement monitoring in this paper. Furthermore, we have also verified the feasibility to extract glacier movement changes and the usability of the PT approach for multi-source imagery. And the evolution of South Inylchek Glacier (SIG) activity will be also discussed simultaneously to some extent. All of these could provide a foundation for the long-term time series analysis of large-scale glacier movement study.

Study area and DataSets
Study area SIG is located in the Tomur Peak region bordering China and Kyrgyzstan, and it is classified as the typical dendritic valley glacier. There is a glacier tongue with an average slope of 2°, which covers an area of approximately 500 km 2 and stretches over 61 km Neelmeijer et al. 2014;Shangguan et al. 2015). The glaciers in the study region begin to melt around in March with the temperature rises, and the annual average temperature is ~ − 7.7 °C (Osmonov et al. 2013). The meltwater of mountain glaciers is the most important freshwater resource for people living downstream because of the Eurasian continental climate with a little rain. Thus, monitoring the SIG movement is of great significance to understand the relationship between glacier activity and climate change in High Mountain Asia (HMA). The detailed image information is shown in Fig. 1, and the glacier outline is described with Randolph Glacier Inventor (RGI) 6.0.

DataSets
In this paper, both optical and SAR imagery acquired during the quasi-same period of the 10-year interval from 2007 to 2008 and 2017 to 2018 was employed in glacier flow estimation respectively, thus the glacier movement calculation would be benefited from multi-source datasets. The optical and SAR data were provided by Landsat-5 and ALOS/PALSAR satellites during 2007-2008 respectively, while the images tracked from 2017 to 2018 were from Landsat-8 and Sentinel-1A. The Near Infrared band of 30 m spatial resolution was employed with better snow penetration performance of Landsat-5 image, while the OLI panchromatic band imaged by Landsat-8 with a 15 m resolution was used (Liu et al. 2020). The relatively high spatial resolution is more conducive to reflecting the details of the glacier surface. Besides, L-band ALOS/ PALSAR and C-band Sentinel-1A SAR images have been exploited for simultaneous comparative analysis of estimated velocity. Compared with the C-band, the L-band SAR image could present stable characteristics for glacier estimation because of its great penetration. Also, both the SRTM DEM and ALOS DEM topographic data were applied to evaluate the elevation changes on the glacier surface, which would be helpful in the analysis of glacier dynamics during the monitoring period. Detailed data information is shown in Table 1. In addition, in order to better reflect the response relationship between glacier movement and climate change in 10 years, the meteorological data in the study area during the study period are shown in Fig. 2. Due to the study area lacks on-site observation data, the temperature and precipitation data of the near weather station (AKSU Station) were selected. At the same time, the comprehensive comparison combined with the dataset of spatially interpolated monthly climate data of WorldClim (Fick and Hijmans 2017), which could assist in the different analysis of long-period glacier movement.

Pixel tracking algorithm
In situ measurement is not appropriate for the investigation of mountain glacier movement in remote regions because it is costly and usually accompanied by danger Wychen et al. 2021). To obtain sub-pixel glacier offset, the pixel tracking algorithm was generally used to estimate the movement of two satellite images on the glacier surface (Huang and Li 2011;Yan et al. 2016;Liu et al. 2017). In this paper, the Normalized Cross-Correlation (NCC) algorithm was utilized in image intensity information for offset computation (Debella-Gilo and Kääb 2011), and the corresponding schematic diagram is shown below (Fig. 3).
In the initial matching process, the template window on the master image was selected firstly, and then the template window slides continuously in the search window of the slave image, which could calculate the normalized correlation coefficient between the template window and the search window until the whole image is traversed. The window with the largest cross-correlation coefficient is the best matching window. Finally, the total displacement of the feature points in two windows is calculated for the primary registration. In order to obtain sub-pixel glacier offset, an oversampling operation is usually undertaken in the matching operation, which could be done by using the two-dimensional Where, (x, y) is the central coordinate of the template window, f is the template window in the master image, g is the search window in the slave image, u and v are the offsets between the two windows in the x-and y-directions, respectively. f is the mean value of the template window in the master image, g is the mean value of the search window in the salve image NCC method uses the intensity information of satellite images to estimate a total offset, including the glacier movement offset, the global-related offset, the topographic-related offset and noise. Based on the fundamental hypothesis that the non-glacial region should be stable without movement (Neelmeijer et al. 2014;Yan et al. 2019), the glacier movement related offset could be revealed by removing irrelevant offsets. NCC could achieve high accurate registration in the glacier region by setting reasonable patch window and resample parameters (Giles et al. 2009;Huang and Li 2011;Nobakht et al. 2014). Besides, we improve the traditional algorithm with GPU to realize efficient parallel operations, which would guarantee the effective monitoring of the glacier movement. The data processing flowchart based on PT is shown in Fig. 4.

Remove irrelevant information
The global related offset was mainly brought by the difference in satellite attitude and position during the imaging process. Generally, the quadratic polynomial fit is capable of covering the global deformation (Yan et al. 2019). To precisely obtain the deformation model, the Random Sample Consensus (RANSAC) algorithm was adopted for a more accurate calculation of the coefficient. During the observation period, we assumed that there was no displacement in most regions of the image, then the transformation model parameters would be efficiently revealed and the orbitrelated offset could be accurately removed.
Besides, the topographic-related offset would appear due to the different imaging angles in mountain regions, and this phenomenon might be serious in SAR images because of its oblique project. To be specific, the azimuth offset is greatly impacted by the intersection angle for orbits and the sensor incidence angle, while the offset along the range direction is proportional to the perpendicular baseline. For SAR imagery, the offset along the azimuth and range directions should be calculated for terrain compensation and an external DEM (Riveros et al. 2013;Yan et al. 2013). And orthorectification would be utilized for optical imagery to remove the terrain displacement. With the assistance of the RPC (Rational Polynomial Coefficient) information of optical images itself and external DEM data, more accurate movement information would be obtained.
After non-glacial movement compensation or removal, the residual and error mismatch points of co-registration might remain. Therefore, the filter operation was needed for reducing the effect of the high-frequency noise. Given both the noise and the edge information belong to the highfrequency signal, the Ranked-order based Adaptive Median Filter (RAMF) was employed for smoothing and de-noising to retain the details (Yan et al. 2019). RAMF works well in glacier movement estimation, as it could not just ensure the accuracy of offset, but also properly filter the noise.

Velocity decomposition
In terms of a geometric difference of projection between SAR and optical data, to make a consistent comparison, all images we used must refer to the same coordinate system before extracting the glacier velocity, which is an essential geometric basis for the comparative analysis of multi-source imagery. To deliver a clearer analysis, the SAR slant range image was decomposed into the plane coordinate system. Given that the terrain of the study area is basically even and smooth gently, we assume that the glacier movement was parallel to the local slope and consistent with the direction in which this maximum slope occurred. The geometric relationship model of SAR image forming is shown in Fig. 5, and the arrow direction is defined as the positive direction in the figure (Ng et al. 2012).
From the geometric relationship of SAR imaging we could get this formula: where D LOS is the surface displacement of the ground point P at the line of sight(LOS), D H and D V are the ground-range (GR) and vertical displacement respectively, and is the sensor incidence angle at the reflection point. In this paper, the (2) D LOS = D H sin − D V cos influence of the glacier deformation in the vertical direction could be reasonably ignored because it has a little contribution to slight displacement (Neelmeijer et al. 2014;Shangguan et al. 2015). D LOS could be regarded as the combination of east-west and north-south displacements on the horizontal plane  where D N is the north-south displacement of the ground point, D E is the east-west deformation, and is the azimuth of the satellite heading vector (positive clockwise from north). To solve the D N and D E , the geometric relationship between the azimuth and the horizontal deformation has been established.
where, D Azi is the ground point deformation along the azimuth direction. With the above formulas, the offset of the SAR image in the east-west and the north-south direction could be obtained.

Results and analysis
The surface movement information of the SIG from 2017 to 2018 was computed based on the Landsat-8 and Sentinel-1A data with the improved PT algorithm as shown in Fig. 6 (3) The Landsat-5 and ALOS/PALSAR data of quasi-same period have been applied in glacier movement estimation of 2007-2008 and corresponding results are shown in Fig. 7. The velocity in the east-west direction proves basically the same despite different data sources. Yet the offset in the north-south direction is less obvious because of the small angle between the direction of the north and the direction in which SAR satellite image, and the inversion results are impacted greatly by noises. Besides, the L-band ALOS/ PASAR data could better reflect the details of the glacier movement because of the excellent penetrability, but the decorrelation of the Landsat-5 image leads to the large differences of velocity in the north-south direction.

Comparative analysis of glacier movement
The total velocity on the glacier surface was obtained by the combination of glacier movement along with different directions (Fig. 8). And the ranges of the color bars were unified (0 cm/d-50 cm/d) for highlighting the features of glacier movement. The velocity changes and elevation differences along the profile line AB (black line in Fig. 1) on the glacier tongue are shown in Fig. 9. A comprehensive analysis of the glacier movement between 2007-2008 and 2017-2018 yields that the velocity distribution looks similar on the profile. Generally, most glacier velocities could achieve 25 to 35 cm/d on the glacier surface. And it decreases gradually from the upper part to the down part. At a distance of ~ 7.5 km from the initial point along profile AB, there is an extreme velocity was observed at the fork with high velocity along the middle line because of the friction from the surrounding mountains. Such velocity distribution shows a satisfying agreement with the result in the previous study Neelmeijer et al. 2014;Nobakht et al. 2014;Yan et al. 2018).
To understand the relationship between the glacier mass variation and the glacier velocity, the "demcoreg" program in ASP was adopted to obtain the elevation change along with the profile on the glacier surface ( Fig. 9) with both SRTM DEM and ALOS DEM (Shean et al. 2016). It indicates that the elevation of the upper and middle reaches has barely changed during around 10 years while the downstream area saw the largest variation. First, the reduction of elevation leads to a decrease of gravitational potential energy (GPE) yet an increase of kinetic energy (KE). Meanwhile, glaciers from the middle and upper reaches of the glacier are constantly transported to the downstream area, and the accumulated stress is released centrally at the fork, this means that the bedrock may be eroded worse at the bottom of the glacier and more likely to move. Finally, the height difference would further increase due to inertia; Second, a large amount of glacial meltwater may flow into Lake Merzbacher during the ablation season, leaving the lake level rising rapidly and much glacier becoming floating after being split off. It may cause the loss of mass at the glacier front, prompting the glacier movement to accelerate for replenishment. A similar monitoring result also could be seen in previous studies Hagg et al. 2018;Li et al. 2018). The accumulation of moraine leads to a slow increase in elevation to the terminus, and the movement velocity decreases gradually when the glacier flows toward the west, which is consistent with the glacier dynamics. The elevation change of the glacier region is relatively small in the study period (Fig. 9), and it could be seen that it is reasonable to ignore the vertical displacement in this paper. Although the research period between the remote imagery and the elevation data is not completely consistent, we found that the elevation curve presents stability in general except for large variation around the turning points (Fig. 9). Therefore, the whole glacier could be generally considered as being in a stable mass balance.
As Figs. 6, 7, 8 show, there is a discontinuity phenomenon in the upstream region of glacier movement fields, which is mainly caused by the saturation associated with the high albedo of clean glacier (Yan et al. 2018). Although the large matching window was selected to perfect the Signal-to-Noise Ratio (SNR), the velocity information might be lost at the edge of the glacier at the same time. In order to analyse the glacier movement better in downstream of the glacier, the flow velocity contours have been shown in Fig. 10. And the ranges of the contour map were unified (5 cm/d to 30 cm/d) for highlighting the features of glacier velocity. Propelled by surface in the middle and upper stream, the glacier branch turning to Lake Merzbacher releases force at the mouth of the lake, causing a stretch at the bottom of the glacier. Here the contours are dense and the glacier velocity changes drastically. There is a large amount of debris on the terminus, which may block the movement of the glacier and be clearly observed in the optical image. In addition, there is a relatively long distance with a gentle slope, causing the glacier velocity to decrease gradually from upstream to downstream. The velocity changes above are consistent with the distribution of moraines as shown in the high-resolution satellite images (Fig. 10).

Accuracy estimation
In this study, the velocity along the same direction respectively from the SAR and optical imagery were compared for the accuracy estimation. And the consistency of the velocity of multi-source images also has been verified. The velocity difference along the profile line in Fig. 1 is shown in Fig. 11. It presents that the velocity fields from multi-source data in different directions are the same in general. Besides, the rationality of the slant-range conversion of SAR data has been verified as the difference of glacier flow rate was ~ ± 2 cm/d in 2017-2018 based on the data of Landsat-8 and Sentinel-1A, which is much smaller than the average velocity. By comparison, the velocity difference could achieve ~ ± 6 cm/d in the north-south direction during 2007-2008. The Near Infrared band of Landsat-5 is insensitive to the glacier variation because of the low resolution, while the L-band ALOS/PALSAR with a full spatial resolution of ~ 5 m could perform better in permeability. Therefore, Fig. 9 Profiles of glacier velocity from different images and corresponding height difference variations there is a large difference in the glacier movement on the glacier surface between the Landsat-5 and ALOS/PALSAR.
Based on the assumption that there is no movement in the glacier-free region in the estimated velocity field, the Root Mean Square Error (RMSE) would be statically computed for evaluating monitoring accuracy (Table 2). And it is roughly the same as the results described by the previous research (Shangguan et al. 2015). The RMSE of the glacierfree region obtained by Landsat-5 images is larger than that of others, taking into account the resolution and other reasons, the experimental result proves valid and acceptable with a high degree of reliability. To ensure the accuracy  of the final monitoring results, the SRTM data has been adopted for the topographic correction. Besides, the satellite images were registered based on the NCC algorithm, which would guarantee the accuracy of co-registration within 0.1 pixels (Ernesto Rodríguez 2006;Mayer et al. 2008).

Discussion
In this study, the SIG velocity has been estimated with both Landsat-8 and Sentinel-1A data in 2017-2018 by the PT algorithm. And the velocity results of different data sources show good consistency in glacier movement. The reliability of the glacier movement estimation approach from multisource was also proved by the comparison of the results in 2007-2008. Obviously, the integration of multi-source imagery could extend the monitoring period and overcome the shortcomings associated with a single data source. However, in the middle part of the glacial trunk, the velocity deviation of the imagery of quasi-synchronization is varied largely. And the corresponding reasons are summarized as follows.
1) Difference in penetration performance. The penetration capabilities of different types of sensors are variable on the glacier surface, and the SAR sensors could achieve deeper than that of the optical sensors. However, there is a little difference in the glacier movement of the quasisame period, which proves that the glacier movements were little affected by the penetration in the study area. Due to the limits of achieved satellite images, the acquisition date of imagery used in glacier movement estimation is not the same, as well as the spatial resolution. Both of them would affect the consistency in glacier velocity distribution of different data sources. Besides, different penetration capability for C-band Sentinel-1A and L-band ALOS/PALSAR would also bring the variation indifference of glacier movements (Strozzi et al. 2008;Yastika et al. 2019;Jefriza et al. 2020). 2) Difference in the research period. The time interval of one year for each image would facilitate the glacier movement calculation and comparison. However, with a different revisit time of satellites and the severely affected optical imagery by weather conditions, the images we could select have more restricted conditions. Meanwhile, the variations in spatial and temporal characteristics happen all-time on the glacier surface, which affects the consistency of the glacier velocity results with a long temporal baseline. 3) Difference in data processing. There is a great difference in terrain correction for SAR and optical imagery due to their different imaging geometry. The optical data was orthorectified with the PRC information and an external DEM for compensating the topographic effect. Meanwhile, with the assistance of DEM and orbit parameters on the PT method, the topographic effect could be accurately compensated in the glacier movement estimation. Besides, the incidence angle of the SAR sensor would change under different topographic regions, and the reliability of the decomposition model in the 3.3 section needs to be verified further.

Conclusions
In this study, the spatio-temporal distributions of glacier movement of SIG were yielded by the improved PT algorithm with optical and SAR imagery during the quasi-same period of the 10-year interval. The results comparison of different datasets with different temporal intervals were given and thoroughly analysed, which would be of great importance to glacier activity and climate change research in Tian Shan. The velocity comparison reveals that there is no significant variation of glacier movement during 2007/2008-2017/2018. Only glacier movement change in the upper and middle parts of the glacier is generally larger than that in the lower part of stream. The rapidly increased movement can be observed at the fork in the downstream, where the maximum elevation change occurs. The observed average glacier velocity ranges from 30 to 35 cm/d, which proves that the movement of glacier surfaces could be accurately estimated by the PT method with both optical and SAR imageries. Furthermore, a combination of imagery acquired by different sensors would be helpful in detecting the glacial movement variation, which would provide an important guarantee for long-term glacier movement analysis. Besides, the RMSE of 0.1 pixels of the average glacier velocity on the glacier-free region verifies the reliability of extracting glacier movement velocity with the multi-source imagery.
With the greatly improved spatial and temporal resolution of SAR and optical imagery, an effective and robust way for both long-term glacier movement analysis and regular glacier resource investigation would be provided. In future research, the meteorological documentation and the accurate DEM data would be also exploited, which could support the analysis of consistent interannual changes in glacier activity along with climate variation in central Asia.