The glacier surface velocity is calculated by the displacement of the characteristics of the glacier surface between different images. For optical image, its displacement can be obtained through the COSI-CORR (Co-Registration of Optically Sensed Images and Correlation), a software package developed based on IDL by California Institute of Technology and then integrated in ENVI software which supports the ENVI version 5 classic (ENVI classic interface). We used the version updated in October 2014 (http://www.tectonics.caltech.edu/slip_history/spot_coseis/download_software.html). The displacement was calculated by COSI-CORR based on the feature matching and mutual correlation of the optical remote sensing image through frequency-domain mutual correlation algorithm.^{14} All phase correlation methods rely on the Fourier transform theory: the relative displacement between a pair of similar images is extracted from the phase difference of their Fourier transform. *i*_{1} and *i*_{2} represent the two images, *x* and *y* are the coordinate values, and *△*_{x} and *△*_{y} represent the displacement of *x* and *y*, respectively, then:

\({i}_{2}\left(x,y\right)={i}_{1}\left(x-{△}_{x},y-{△}_{y}\right)\)**(1)**

where *I*_{1} and *I*_{2} represent the Fourier transform obtained through the Fourier transform theorem, which can be expressed as:

\({I}_{2}\left({\mathrm{\omega }}_{x},{\omega }_{y}\right)={I}_{1}\left({\omega }_{x},{\omega }_{y}\right){e}^{-j\left({\omega }_{x}{△}_{x}+{\omega }_{y}{△}_{y}\right)}\)**(2)**

where *ω*_{x} and *ω*_{y } are the frequency variables in column and row, and *e* is the natural constant, and *j *is the imaginary number. Therefore, the normalized cross-spectrum (*C*) of the images *i*_{1} and *i*_{2 } can be expressed as:

\({C}_{{i}_{1}{i}_{2} }\left({\omega }_{x},{\omega }_{y}\right)=\frac{{I}_{1}\left({\omega }_{x},{\omega }_{y}\right){I}_{2}^{*}\left({\omega }_{x},{\omega }_{y}\right)}{\left|{I}_{1}\left({\omega }_{x},{\omega }_{y}\right){I}_{2}^{*}\left({\omega }_{x},{\omega }_{y}\right)\right|}={e}^{j\left({\omega }_{x}{△}_{x}+{\omega }_{y}{△}_{y}\right)} \)**(3)**

where * denotes the complex conjugate. The Dirac delta function (\(\delta \left(x+{△}_{x},y+{△}_{y}\right)\)) is obtained by inverse Fourier transform (\({\mathcal{F}}^{-1}\)) of the normalized cross-spectrum:

\({\mathcal{F}}^{-1}\left\{{e}^{j\left({\omega }_{x}{△}_{x}+{\omega }_{y}{△}_{y}\right)}\right\}=\delta \left(x+{△}_{x},y+{△}_{y}\right)\)**(4)**

The relative displacement of images can then alternatively be estimated from the coordinates of the correlation peak. In case of subpixel displacements, this peak is not a Dirac delta function anymore, but a down-sampled version of a Dirichlet kernel. Therefore, the phase correlation method of COSI-CORR can be followed by two steps: (1) the relative migration of the images reconstructed by calculating the linear phase of their cross-power spectrum; (2) the relative displacement of the images determined by the position of the strict correlation peak value. In order to calculate the relative displacement, it is necessary to apply ortho-rectification to one of the images, and then use the orthophoto as the reference image to coregister another image. The Landsat data used in this study were obtained from USGS, and all the data were ortho-rectified using the ground control point of GLS2005 by USGS. Guo et al.^{15} verified the correction accuracy, and the results showed that the Landsat image provided by USGS after ortho-rectification had a high correction accuracy: most of the correction precision was around half pixel, and some of the images even reached the accuracy of 1/6 ~ 1/10 pixel. Therefore, it was considered that the accuracy of the ortho-rectified Landsat images satisfied the glacier surface velocity analysis. The Landsat OLI band 8 (panchromatic band, with a spatial resolution of 15 m) was calculated by the COSI-CORR, where frequency domain algorithm was used to calculate the correlation coefficient, with the reference window set to 128 and the search window set to 32. The resulting displacement data consisted of three images that recorded east-west displacement (EW), north-south displacement (NS) and signal-to-noise ratio (SNR), respectively. We removed the results of SNR < 0.9 and cloud or shadow cover, and obtained the glacier surface velocity data.