Ice, Snow and Environment Over High Asia Zone II Versions EN1 Vol 3 (4) 2018
Dataset of glacier surface motion along KKH during 2015 – 2017
>>
： 2018 - 03 - 19
： 2019 - 03 - 27
： 2018 - 11 - 09
617 5 0
Abstract & Keywords
Abstract: As an international highway, the Karakoram Highway (KKH) is of great importance to both China and Pakistan. Passing through one of the most glacier-concentrated regions—the Karakoram Mountains, KKH constantly faces threats from glacier motions along the way. This study obtains Sentinel-1A SAR data of the study area, and uses the feature-tracking method to calculate the field of glacier surface displacement along the KKH during October 2014 and March 2018. Results indicated that while glacier surface velocities in the northern and central parts of the region were relatively stable, some large-scale glaciers in the southern part had higher velocities up to 4 m/d, particularly along the section near Batura Sar. This dataset can be used in research of glacial mass balance, hydrological resource dynamics and climate changes. In addition, the high-precision monitoring results of the displacement field are of great significance for glacier disaster warnings to KKH and the residents around.
Keywords: KKH; feature-tracking; glacier surface motion; Sentinel-1A
Dataset Profile
 Chinese title 2015–2017年中巴公路沿线山地冰川表面运动数据集 English title Dataset of glacier surface motion along KKH during 2015 – 2017 Data corresponding authors Jiang Zongli (jiangzongli@hnust.edu.cn) Data authors Wang Lei, Jiang Zongli, Liu Shiyin, Zhang Zhen Time range October 31 2014 to March 8 2018 (with a 12/24-day interval) Geographical scope 35°50′53″N – 39°00′42″N, 74°17′58″E – 75°53′01″E Spatial resolution 100m × 100m Data volume 9.51 GB Data format *.tif (GeoTIFF, 32 bit float) Data service system Sources of funding National Natural Science Foundation of China (No. 41471067, 41474014, 41701087, 41761144075); Special Program of the Ministry of Science and Technology (2013FY111400) Dataset composition The dataset consists of two data files: “Direction” stores raster data of glacier surface motion directions, and “Velocity” stores data of glacier surface motion speeds.
1.   Introduction
Karakoram Highway (KKH) is an international highway significant for economic and cultural exchange among China, Pakistan and other countries along the way. It is not only an important witness of the friendship between China and Pakistan, not also provides a vital channel for communications between western China and central, western and south Asia. The road passes through the Karakoram mountains which host the largest number and scale of glaciers.1 KKH has a typical continental climate featuring small precipitation and great temperature variation in summer and winter. In comparison, the glaciers along KKH were located at a high altitude featuring typical mountain climate. Effected by weakened Indian Ocean Monsoon and strengthened Westerly, certain high altitude glacier areas had an annual precipitation of even more than 2000mm·a-1.2
There are more than 2500 glaciers assembled within 50 km of KKH, including large-scale high-mass-balance glaciers like Hispar Glacier developed on K2, Batura Glacier and Pasu Glacier developed on Batura Sar.3 The glacier disasters imperil the highway all the year around. In 2008, Ghulkin glacier’s accelerated advancement resulted in 4 debris flows triggered by glacier lake outbursts, causing serious damage to the highway.4 In 2015, Amaury Dehecq used optical RS images to obtain the glacier surface velocity in the Karakorum area, finding that most of the glaciers kept stable motions with exceptions of large-scale glaciers that had a faster velocity.5-7 Many researches showed abnormal motions of the glacier surface, though variations occurred between findings of different periods.8,9 These conclusions proved great complexity of the glacier surface motion in this area.

Fig.1   Location of KKH and the glaciers nearby
With its continuous earth observation capability over large-scale areas, remote sensing has gradually replaced the traditional methods of surveying and mapping like GNSS (Global Navigation Satellite System) in the monitoring of extreme environment like glacial area. For decades now, SAR (Synthetic Aperture Radar) enables multiband, multi-polarized and multi-mode continuous data acquisition free from climatic or atmospheric interference,10,11 which can be used to observe large-scale glacier surface motion.
With longer intervals, SAR image pairs have problems like relatively serious temporal decorrelation, while feature-tracking method, by using pixel-by-pixel matching to search characteristic points, can minimize the negative impact of temporal decorrelation. In recent years, feature-tracking method has been used to get the surface motion data of glaciers worldwide. Particularly in 2002, Strozzi et al. used it to extract the glacier surface motion velocity, which proved its feasibility and reliability in glacier surface motion analysis.12,13 In this study, we obtained 63 phase Sentinel-1A data along the KKH during October 2014 and March 2018, used the feature-tracking method to analyze the glacier surface motion, and obtained 58 pairs of successive glacier surface motion data.
2.   Data collection and processing
2.1   Data sources
Developed by ESA (European Space Agency), Sentinel-1A is a full-time full-coverage earth observation imaging system equipped with synthetic aperture radar. As shown in Figure 1, the study area has a strip-like terrain and needs to be covered by three Sentinel-1A images. Our study selected data under IW mode and VV polarization mode from three consecutive positions in Path 27, consisting of 189 images obtained during October 2014 and March 2018 (63 images for each position). Table 1 shows some details of the data.
Table 1   Sentinel-1A data used in this study
 Position Time pair Image count Time interval (day) Polarization mode Observation mode Path 27 frame 121 2014-10-06～2016-12-31 30 24 VV IW 2017-01-24～2018-03-08 33 12 Path 27 frame 117 2014-10-06～2016-12-31 30 24 2017-01-24～2018-03-08 33 12 Path 27 frame 112 2014-10-06～2016-12-31 30 24 2017-01-24～2018-03-08 33 12
Notes: VV, Vertical like-Polarized; IW, Interferometric Wide swath observation mode; the time interval of image pair during 20150710-20150827 is 48 days, and that during 20150827-20151002 is 36 days.
Additionally, we obtained 30 m resolution SRTM (Shuttle Radar Topography Mission) DEM from NASA (National Aeronautics and Space Administration) and USGS (United States Geological Survey), based on which terrain-based precise geometric correction was performed on SAR data and calculation to extract the altitude and slope of research area.
2.2   Methods and procedure
Glacier is known to be fluid. Particularly, the surface motion of a mountain glacier can reach a speed of hundreds of meters a year during its surging period. It means that the SAR image pairs of a glacial area are of low correlation. Therefore, feature-tracking method was used in this study to extract the shifting of feature point between SAR image pairs as it was less affected by the low correlation of the image pairs.
Feature-tracking method for SAR image is theoretically similar to that for optical image which conducts global feature searching based on a given search window for obtaining the feature area and feature point matching pairs, and for eliminating the feature point matching pairs overly distant from the given offset threshold. Afterwards, the feature point’s component offset on ranges and azimuth is calculated by using the offset-tracking method and the maximum correlation coefficient of image texture. Then the 2D offset vectors of SAR images’ feature point pairs are composed by using its component offset vectors,14,15 which denote the offset direction of the SAR image’s feature point and the offset value on the direction. The registration offsets contain Roffset (motion offsets of the same ground feature point at different times, orbital offsets Rorbit (orbit deviation offsets caused by imprecise repeated orbit when the same satellite captures the image of the same object), and Rorbit (offsets caused by rugged terrain). The registration offsets $${R}_{offset}$$ can be expressed by:
${R}_{offset}={R}_{motion}+{R}_{orbit}+{R}_{dem}$
The offsets caused by rugged terrain are small enough to be ignored in the image pairs with a shorter baseline or flat terrain. So the registration offsets can be expressed by:
${R}_{offset}={R}_{motion}+{R}_{orbit}$
It can be considered that in the non-glacier area of this study, there are no offsets on image pairs with shorter time baselines. In order to extract the offsets of glacier feature point, we need to remove the orbital offsets. So we oversample the matching windows, and fit the orbital offsets using the least square function:
$R={a}_{0}+{a}_{1}x+{a}_{2}y+{a}_{3}xy$
where $${a}_{0}，{a}_{1}，{a}_{2}，{a}_{3}$$ are indeterminate coefficients, $$x\mathrm{ }$$expresses the value of range in the SAR data coordinate system, $$y$$ expresses the value of azimuth in the SAR data coordinate system. Emulate the function by fitting the peak values determined by interpolation.
Gamma and ArcGIS were used to process and analyze the data. Figure 2 shows the specific procedure. Before the calculation, SRTM data were used to perform registration on SAR images to compensate the errors caused by rugged terrain. In order to improve accuracy, non-glacial areas were first extracted by visual interpretation based on SAR images, and then a search window of 128×256 pixels was used to match the feature point and feature area. Finally, the feature-tracking method was used to calculate the offsets of theses points and areas.

Fig.2   Schematic diagram illustrating the glacier surface motion velocity
The procedure consists of four steps: preprocessing, images registration, offsets estimation and coordinate transformation. In preprocessing, the Sentinel-1A level 1 data and SRTM DEM were converted into a Gamma-processable format, and meanwhile, some table and parameter files were prepared for the next step. In images registration, SAR registration theory was used to extract the images’ feature point offsets caused by the two-dimensional deformation of earth surface, on which topographical correlation was performed to refine offset information assisted by SRTM data. Then the main-images and slave-images were roughly registered based on the global displacement of orbit parameters, and the initial offset with an accuracy of 0.5 pixels was obtained. According to the initial offset and the pixel-level registration principle, a certain size was set for the search window and a certain length was set for the search step, based on which the precise offset of similar surfaces was obtained by calculating the coherence. The matching search-window of 128×256 pixels was used in this study.
The estimation of offsets was based on the feature matching results of image pairs obtained in two time periods. The standard deviation was taken as the accuracy of registration. Surface displacement was obtained by subtracting the fitted orbit deviation from the registration offsets. Since the sensor of Sentinel-1A had a certain incidence angle, there were some shadows or overlaps which led to invalid data. This was more conspicuous in the terrain of large fluctuations like glacier-covered area. Finally, slope data were introduced for calculating 3D movement vectors of the glacial area.
3.   Sample description
The dataset of glacier surface motion along KKH during 2015 – 2017 mainly contains two parts of data. One is for the direction of the glacier surface motion in the study area, named “YYYYMMDD_YYYYMMDD_direction.tif”. It is expressed by radian and stored in the folder “Direction”, recording the movement direction of the glacial feature points during two time periods. The other is named “YYYYMMDD_YYYYMMDD_velocity.tif”, which records the motion velocity of glacier surface feature points in the research area. It is stored in the folder “Velocity”. The data are in the unit of m/d. Sample data are shown in Figure 3. Figure 3(a) shows the annual average velocities of glacier surface motion around Batura Sar, while Figure 3(b) shows the annual average velocities of glacier surface motion around Mount Kongur.

Fig.3   Glacier surface velocities along the KKH
4.   Quality control and assessment
The Sentinel-1A data in IW mode this dataset used had a resolution of 5m×20m, which might ineluctably cause data errors. In addition, errors also occurred in the data calculation process like images registration. The errors were mainly from the offset estimation of the image pairs, approximately equal to the sum of the errors of the range and azimuth directions.16 The error of the image registration algorithm, provided by the Swiss GAMMA radar data processing platform, was controlled within 0.01 pixels, and the corresponding range and azimuth directions of Sentinel-1A were about 0.018 m and 0.009 m respectively. The errors of this magnitude had limited influence on the glacier motion monitoring.
Table 2   Statistics showing errors in the non-glacier area (m/d)
 Data pairs Variance Data pairs Variance Data pairs Variance 20141031_20141124 0.0540 20160809_20160902 0.0048 20170711_20170723 0.0835 20141218_20150111 0.0183 20160902_20160926 0.0042 20170723_20170804 0.1271 20150111_20150204 0.0414 20160926_20161020 0.0017 20170804_20170816 0.0744 20150312_20150405 0.0662 20161020_20161113 0.0012 20170816_20170828 0.0632 20150405_20150429 0.0917 20161113_20161207 0.0017 20170828_20170909 0.0606 20150429_20150523 0.0785 20161207_20161231 0.0013 20170909_20170921 0.0660 20150523_20150616 0.0680 20161231_20170124 0.0029 20170921_20171003 0.0531 20150616_20150710 0.0833 20170124_20170205 0.0202 20171015_20171027 0.0698 20150710_20150827 0.0300 20170205_20170217 0.0079 20171027_20171108 0.0630 20150827_20151002 0.0419 20170217_20170301 0.0101 20171108_20171120 0.0604 20151002_20151119 0.0406 20170301_20170313 0.0175 20171120_20171202 0.0864 20151213_20160106 0.0466 20170313_20170325 0.0167 20171202_20171214 0.1005 20160106_20160130 0.0022 20170325_20170406 0.0543 20171214_20171226 0.0622 20160130_20160223 0.0019 20170406_20170418 0.1159 20171226_20180107 0.0610 20160223_20160318 0.0032 20170418_20170430 0.0829 20180107_20180119 0.0558 20160318_20160411 0.0176 20170430_20170512 0.0406 20180119_20180212 0.0216 20160411_20160505 0.0346 20170512_20170524 0.0443 20180212_20180224 0.0676 20160505_20160529 0.0119 20170524_20170605 0.0206 20180224_20180308 0.0786 20160529_20160716 0.0041 20170605_20170629 0.0556 Average 0.0445 20160716_20160809 0.0015 20170629_20170711 0.0816
In addition to the errors caused by the GAMMA software mentioned above, other error sources of this dataset include image registration error, terrain-induced ground migration, geocoding-induced error, orbital error, and so forth. In order to quantify the errors, we assumed that the non-glacier area remains static within 12 or 24 days and thus the resulting offsets in the non-glacier area include all of the above error sources. We calculated statistics in the non-glacier area, which showed a holistic variance of 0.0445 m/d. Details are shown in Table 2. Figure 4 shows the histogram of the holistic sampling errors. It can be seen that the accuracy of the resulting data suffices a spatial resolution of 100m * 100m. As mentioned above, the sampling error analysis assumed all non-glacier areas as static, whereas in fact, small-scale changes took place in the areas, like small water area changes and small surface changes resulting from snow melting at piedmont in spring and summer. The errors of the sampling statistics in the non-glacier area hence inevitably represent the maximal error of the resulted data.
For high mountains in the study area, radar sensor sometimes failed in observing the overlapping and shadowing areas, which led to less credible calculation results. In these circumstances, we replaced these results with the data extracted by using viewshed analysis. Furthermore, because of misregistration when the offset distance of feature point matching pairs exceeded the set offset threshold, a small scale of data holes appeared in the glacier accumulation area in late spring and early summer. The reasons could be attributed to drastic surface temperature changes, strong wind, steep slope and avalanche apt to take place in later spring and earlier summer.

Fig.4   Errors of glacier motion velocity in the non-glacier area
5.   Data value and application
Using the C band synthetic aperture radar sensor, Sentinel-1A was able to minimize the impact of atmospheric interference. Additionally, with the shorter time interval of its images, the glacier motion data achieved a higher precision and time resolution. For a better calculation precision, SRTM DEM was used for offsetting the topographic migration in the procedure of SAR images registration. Glacier surface motion data are fundamental for studying glacier mass balance, glacier movement characteristics and early warning of glacier disasters. When used in combination with climatic and hydrological data, the glacier surface motion data can help comprehensively analyze the characteristics and trends of regional geographic environment changes. The results of glacial mass balance analysis are of great significance to the study of the distribution and change of regional and global water resources.
6.   Usage notes and recommendations
The dataset of glacier surficial motion along KKH during 2015 – 2017 provides whole data download service through CERN Open Data Portal. This dataset provides 32 bit floating point grid data in GeoTIFF format. Flow velocity data are expressed as pixel value in the unit of m/d, and the flow direction data are stored as pixel value in radian. This dataset can be read by GIS software like ArcGIS and SuperMap, or processed and analyzed in batches using data analysis software like Matlab.
Acknowledgments
The authors would like to thank ESA (European Space Agency) for supplying Sentinel 1A data, CGIAR-CSI (Consultative Group for International Agricultural Research Consortium for Spatial Information) for supplying 30m resolution SRTM SEM data, GLIMS (Global Land Ice Measurements from Space) and Scientific Data Center, Institute of Cold and Arid Regions, Chinese Academy of Sciences for supplying glacier inventory dataset.
1.
Zhang X, Chen JM & Wang WY. Recent variations of the Batura Glacier in the Karakorum Mountains. Journal of Glaciology and Cryopedology 18(1996): 13.
2.
Shi Y & Zhang X. Historical variations in the advance and retreat of the Batura Glacier in the Karakoram Mountains. ACTA Geographica Sinica 33(1978): 16.
3.
Shi Y, Ding D & Fu L. The numerical forecasting model of the terminal changes in the Batura Glacier. Journal of Glaciology and Cryopedology 4(1982): 10.
4.
Liao L, Zhu Y, Yang Z et al. Advance and Retreat Fluctuation of the Ghulkin Glacier Along the Karakoram Highway over Hundred years. Journal of Glaciology and Geocryology 35(2013): 9.
5.
Akhtar M, Ahmad N & Booij MJ. The impact of climate change on the water resources of Hindukush–Karakorum–Himalaya region under different glacier coverage scenarios. Journal of Hydrology 355(2008): 148 – 163.
6.
Ashraf A, Naz R & Iqbal MB. Altitudinal dynamics of glacial lakes under changing climate in the Hindu Kush, Karakoram, and Himalaya ranges. Geomorphology 283(2017): 72 – 79.
7.
Muhammad S & Tian L. Changes in the ablation zones of glaciers in the western Himalaya and the Karakoram between 1972 and 2015. Remote Sensing of Environment 187(2016): 505 – 512.
8.
Khromova TE, Osipova GB, Tsvetkov DG et al. Changes in glacier extent in the eastern Pamir, Central Asia, determined from historical data and ASTER imagery. Remote Sensing of Environment 102(2006): 24 – 32.
9.
Dehecq A, Gourmelen N & Trouve E. Deriving large-scale glacier velocities from a complete satellite archive: Application to the Pamir–Karakoram–Himalaya. Remote Sensing of Environment 162(2015): 55 – 66.
10.
Yang K, Yang J & Jiang B. Sentinel-1 satellite overview. Urban Geotechnical Investigation & Surveying (2015): 24 – 27.
11.
Zhang Z, Liu S, Wei J et al. Monitoring recent surging of the Karayaylak Glacier in Pamir by remote sensing. Journal of Glaciology and Geocryology 38(2016): 11 – 20.
12.
Shangguan D, Liu S, Ding Y et al. Characterizing the May 2015 Karayaylak Glacier surge in the eastern Pamir Plateau using remote sensing. Journal of Glaciology 62(2016): 944 – 953.
13.
Li R, Ma X, Cheng Y et al. Ice velocity measurement in east Antarctica from 1960s to 1980s based on Argon and Landsat imagery. ISPRS - International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences (2017): 1535 – 1539.
14.
Li J, Li Z, Wang C et al. Using SAR offset-tracking approach to estimate surface motion of the South Inylchek Glacier in Tianshan. Chinese Journal of Geophysics 56(2013): 11.
15.
Jiang Z, Liu S, Xu J et al. Using feature-tracking of ALOS PALSAR images to acquire the Yengisogat glacier surface velocities. Journal of Glaciology and Geocryology 33(2011): 512 – 518.
16.
Qureshi MA, Yi C, Xu X et al. Glacier status during the period 1973–2014 in the Hunza Basin, Western Karakoram. Quaternary International 444(2017): 125 – 136.
Data citation
1. Wang L, Jiang Z, Liu S et al. Dataset of glacier surface motion along KKH during 2015 – 2017. Science Data Bank. DOI: 10.11922/sciencedb.628 (2018).
Article and author information
Wang L, Jiang Z, Liu S et al. Dataset of glacier surface motion along KKH during 2015 – 2017. China Scientific Data 4(2018). DOI: 10.11922/csdata.2018.0006.zh
Wang Lei
SAR data processing.
postgraduate student; research area: glacier surface changes.
Jiang Zongli
SAR data processing schemes.
jiangzongli@hnust.edu.cn
PhD, Associate Professor; research area: frozen zone change based on remote sensing technology.
Liu Shiyin
SAR data processing schemes.
Professor, doctoral supervisor; research area: frozen zone changes.
Zhang Zhen
processing flowchart design.
PhD; research area: frozen zone change based on remote sensing technology.
National Natural Science Foundation of China (No. 41471067, 41474014, 41701087, 41761144075); Special Program of the Ministry of Science and Technology (2013FY111400)
Publication records
Published: Nov. 9, 2018 （ VersionsEN1
Released: March 26, 2018 （ VersionsZH2
Published: Nov. 9, 2018 （ VersionsZH3
References

csdata