Ice Velocity Variations of the Polar Record Glacier ( East Antarctica ) Using a Rotation-Invariant Feature-Tracking Approach

In this study, the ice velocity changes from 2004 to 2015 of the Polar Record Glacier (PRG) in East Antarctica were investigated based on a feature-tracking method using Landsat-7 enhanced thematic mapper plus (ETM+) and Landsat-8 operational land imager (OLI) images. The flow field of the PRG curves make it difficult to generate ice velocities in some areas using the traditional normalized cross-correlation (NCC)-based feature-tracking method. Therefore, a rotation-invariant parameter from scale-invariant feature transform (SIFT) is introduced to build a novel rotation-invariant feature-tracking approach. The validation was performed based on multi-source images and the making earth system data records for use in research environments (MEaSUREs) interferometric synthetic aperture radar (InSAR)-based Antarctica ice velocity map data set. The results indicate that the proposed method is able to measure the ice velocity in more areas and performs as well as the traditional NCC-based feature-tracking method. The sequential ice velocities obtained present the variations in the PRG during this period. Although the maximum ice velocity of the frontal margin of the PRG and the frontal iceberg reached about 900 m/a and 1000 m/a, respectively, the trend from 2004 to 2015 showed no significant change. Under the interaction of the Polar Times Glacier and the Polarforschung Glacier, both the direction and the displacement of the PRG were influenced. This impact also led to higher velocities in the western areas of the PRG than in the eastern areas. In addition, elevation changes and frontal iceberg calving also impacted the ice velocity of the PRG.


Introduction
The contribution of mass loss from the East Antarctic ice sheet over the past 20 years remains unclear [1].Ice flows and calving cause the rapid advance and retreat of marine-terminating glaciers, therefore outlet glaciers often respond more rapidly to the local effects of global climate change that inland glaciers [2].Therefore, for a better understanding of the ongoing changes of mass balance in East Antarctica, studies of outlet glaciers and their changes in ice velocity need to be conducted.
A number of studies of outlet glaciers in East Antarctica have been conducted.Stearns et al. [3] provided direct evidence for the acceleration of the ice velocity on the Byrd Glacier in response to two subglacial lake drainage events.Golledge and Levy [4] suggested that the ice velocity of the East Antarctic Ice Sheet outlet glaciers has responded dynamically in response to environmental changes.Miles et al. [5] analyzed the calving activity of six outlet glaciers in the Porpoise Bay region (East Antarctica) between November 2002 and March 2012, and investigated the drivers behind the observed calving dynamics.In addition, observed variations in the ice velocities of different outlet glaciers and ice shelves in East Antarctica have been reported, including the Amery Ice Shelf [6], the Mertz Glacier Tongue [7], the Langhovde Glacier [8], the Shirase Glacier [9], and the Stancomb-Wills Ice Tongue [10].The Polar Record Glacier (PRG) is the largest outlet glacier along the Ingrid Christensen Coast in East Antarctica.Zhou et al. [11] presented a study of the seasonal and interannual ice velocity changes at the PRG in 1996PRG in , 2006PRG in , 2007PRG in , and 2008.However, there have been few studies of the long-term and annual ice velocity variations of the PRG.
Remotely-sensed imagery can provide detailed and timely data for Earth observation across both time and space, which enhances the ability to map and monitor glacier flow on a nearly global scale.Optical-satellite-image-based ice-velocity measurement using feature tracking is a well-established method.Feature tracking involves tracking identifiable features between pairs of optical satellite images using an image-matching algorithm such as normalized cross-correlation (NCC) [12][13][14][15][16][17], cross-correlation operated in the frequency domain on orientation images (CCF-O) [18][19][20], and co-registration of optically sensed images and correlation (COSI-Corr) [21].The results of a comparison [22] indicate that COSI-Corr is the most robust matching method.However, COSI-Corr has problems on narrow outlet glaciers [22].Although NCC performs poorly in areas of poor visual contrast and areas with changing snow conditions, it performs better than the other methods on narrow outlet glaciers such as the PRG [22].
Therefore, in this research, ice velocity changes at the PRG from 2004 to 2015 were investigated using the NCC-based feature-tracking method with Landsat series images.However, the curved PRG makes it difficult to generate estimates in some areas using the traditional NCC-based feature-tracking method.Scale-invariant feature transform (SIFT) [23], which transforms image data into scale-and orientation-invariant coordinates relative to local features, has proven its efficiency for various kinds of applications in optical remote sensing [24].Mikolajczyk and Schmid's research also presents the robustness and the distinctive character of the SIFT [25].In this method, each keypoint is assigned a consistent orientation based on local features and described relative to this orientation, by which invariance to the image rotation can be achieved.Thus, the rotation-invariant parameter from SIFT, which is used to extract the consistent orientation, is introduced into the traditional NCC-based feature-tracking method to make it suitable for feature tracking of the PRG.
We briefly introduce the PRG and the data sources in Section 2. In Section 3, the rotation-invariant feature-tracking method is introduced, and its accuracy is validated.The ice velocity variations of the PRG are analyzed and discussed in Sections 4 and 5. Section 6 concludes the research.

Study Area
The PRG is located in the Prydz Bay area (Figure 1) and is the largest outlet glacier along the Ingrid Christensen Coast.The widest part of the glacier reaches about 17 km, and the narrowest place is about 8.5 km.The Polar Times Glacier, the Polarforschung Glacier, and the Amery Ice Shelf lie to the west, and Zhongshan Station lies to the east.In front of the PRG, there is an iceberg that broke away from the glacier between 1973 and 1989 [11].

Sequential Landsat Series Images and Preprocessing
Landsat images are easily-acquired, high-resolution optical satellite data for polar regions.The main advantages of Landsat images are that one image covers a large area.Moreover, data from this sensor extend back to 1972, making it possible to conduct a long-term study.Fifteen Landsat-7 enhanced thematic mapper plus (ETM+) and Landsat-8 operational land imager (OLI) images were used for the PRG monitoring.These images were provided by the United States Geological Survey and chosen at approximately annual intervals (Table 1).The spatial resolution of the multispectral bands in each image is 30 m.Since the scan line corrector (SLC) in the ETM+ instrument failed on 31 May 2003, stripes have existed in the Landsat-7 ETM+ images.Therefore, destriping was performed using a local histogram matching method developed by the National Aeronautics and Space Administration (NASA) [26] before the ice velocity derivation.

Sequential Landsat Series Images and Preprocessing
Landsat images are easily-acquired, high-resolution optical satellite data for polar regions.The main advantages of Landsat images are that one image covers a large area.Moreover, data from this sensor extend back to 1972, making it possible to conduct a long-term study.Fifteen Landsat-7 enhanced thematic mapper plus (ETM+) and Landsat-8 operational land imager (OLI) images were used for the PRG monitoring.These images were provided by the United States Geological Survey and chosen at approximately annual intervals (Table 1).The spatial resolution of the multispectral bands in each image is 30 m.Since the scan line corrector (SLC) in the ETM+ instrument failed on 31 May 2003, stripes have existed in the Landsat-7 ETM+ images.Therefore, destriping was performed using a local histogram matching method developed by the National Aeronautics and Space Administration (NASA) [26] before the ice velocity derivation.
Accurate image registration is required in the time series analysis of velocity fields.In addition, radiance calibration, atmospheric correction, and illumination correction are also required in the image preprocessing.In this research, ENVI 4.6 software was used for the image preprocessing.

Validation Data
The making earth system data records for use in research environments (MEaSUREs) interferometric synthetic aperture radar (InSAR)-based Antarctica ice velocity map data set [27] was used as the validation data.This data set provides two resolutions of ice velocities for Antarctica: 450 m and 900 m.The resolution used in this research was 450 m.The error for this data set is from 1 m/a to greater than 16 m/a [27].The majority of the data set was acquired from multi-source satellite images obtained from 2007 to 2009.

Rotation-Invariant Feature Tracking
Tracking features between pairs of images is a robust way to define regions of movement such as the ice flow in Antarctica [13,15].This method can be used with both optical satellite images and also synthetic aperture radar satellite images [7,8,28], and has been proven to be effective for most glaciers.
Figure 2a,b present the ice velocity maps for the Polar Times Glacier, which rotates with a small angle after it calves.Figure 2c presents the ice velocity map which was derived using the traditional NCC-based feature-tracking method.Figure 2c indicates that this method is unable to generate estimates in some areas on the front of the outlet glacier.The reason for this can be explained using Figure 3.
In Figure 3, we refer to the red pixel in the slave image as the matching point which corresponds to the keypoint in the master image.M is the reference chip in the master image.S and S are match chips in the slave image.In a non-curved glacier, M will directly flow to S over the time interval.Once the glacier curves, M will flow to S and rotate at an orientation angle θ compared with S. The traditional NCC-based feature-tracking method can obtain the displacements in both the x-axis and y-axis [28], but neglects the rotated variation coming from the angle θ.This leads to the result that the ice velocities in some areas could not be generated.Therefore, the traditional NCC-based feature-tracking method is not suitable for a curved outlet glacier.
Accurate image registration is required in the time series analysis of velocity fields.In addition, radiance calibration, atmospheric correction, and illumination correction are also required in the image preprocessing.In this research, ENVI 4.6 software was used for the image preprocessing.

Validation Data
The making earth system data records for use in research environments (MEaSUREs) interferometric synthetic aperture radar (InSAR)-based Antarctica ice velocity map data set [27] was used as the validation data.This data set provides two resolutions of ice velocities for Antarctica: 450 m and 900 m.The resolution used in this research was 450 m.The error for this data set is from 1 m/a to greater than 16 m/a [27].The majority of the data set was acquired from multi-source satellite images obtained from 2007 to 2009.

Rotation-Invariant Feature Tracking
Tracking features between pairs of images is a robust way to define regions of movement such as the ice flow in Antarctica [13,15].This method can be used with both optical satellite images and also synthetic aperture radar satellite images [7,8,28], and has been proven to be effective for most glaciers.
Figure 2a,b present the ice velocity maps for the Polar Times Glacier, which rotates with a small angle after it calves.Figure 2c presents the ice velocity map which was derived using the traditional NCC-based feature-tracking method.Figure 2c indicates that this method is unable to generate estimates in some areas on the front of the outlet glacier.The reason for this can be explained using Figure 3.
In Figure 3, we refer to the red pixel in the slave image as the matching point which corresponds to the keypoint in the master image.M is the reference chip in the master image.S and Sʹ are match chips in the slave image.In a non-curved glacier, M will directly flow to S over the time interval.Once the glacier curves, M will flow to Sʹ and rotate at an orientation angle compared with S. The traditional NCC-based feature-tracking method can obtain the displacements in both the x-axis and y-axis [28], but neglects the rotated variation coming from the angle .This leads to the result that the ice velocities in some areas could not be generated.Therefore, the traditional NCC-based featuretracking method is not suitable for a curved outlet glacier.To address this issue, before feature tracking, the master image should first be rotated at an orientation angle .As a result, by keeping the same orientation as the slave image, the proposed method achieves invariance to image rotation.The process flow for rotation-invariant featuretracking is presented in Figure 4.The fundamental principle of the feature tracking method is that persistent surface features common to sequential images can be tracked and identified.The corresponding spatial displacements of these features indicate the ice velocity.The procedures of the rotation-invariant feature-tracking method are described in detail as follows: 1. Preprocessing.The master image and slave image are strictly georeferenced and orthorectified.
By performing this step, the relief displacement and geometric distortion can be removed.2. Keypoint localization.In master image, keypoints were chosen at regular intervals.Since the resolution of the MEaSUREs InSAR-Based Antarctica Ice Velocity Map used in this research is 450 m, the corresponding resolution of the Landsat series images is 30 m.Therefore, the interval is 13 pixels.3. Orientation assignment.The rotation-invariant parameter coming from SIFT is used to extract orientation angle .Firstly, the gradient magnitude ( , ) and orientation ( , ) are the bases for the calculation of orientation .They can be calculated utilizing the features from the keypoint and its neighborhoods, as follows [23]: To address this issue, before feature tracking, the master image should first be rotated at an orientation angle θ.As a result, by keeping the same orientation as the slave image, the proposed method achieves invariance to image rotation.The process flow for rotation-invariant feature-tracking is presented in Figure 4. To address this issue, before feature tracking, the master image should first be rotated at an orientation angle .As a result, by keeping the same orientation as the slave image, the proposed method achieves invariance to image rotation.The process flow for rotation-invariant featuretracking is presented in Figure 4.The fundamental principle of the feature tracking method is that persistent surface features common to sequential images can be tracked and identified.The corresponding spatial displacements of these features indicate the ice velocity.The procedures of the rotation-invariant feature-tracking method are described in detail as follows: 1. Preprocessing.The master image and slave image are strictly georeferenced and orthorectified.
By performing this step, the relief displacement and geometric distortion can be removed.2. Keypoint localization.In master image, keypoints were chosen at regular intervals.Since the resolution of the MEaSUREs InSAR-Based Antarctica Ice Velocity Map used in this research is 450 m, the corresponding resolution of the Landsat series images is 30 m.Therefore, the interval is 13 pixels.3. Orientation assignment.The rotation-invariant parameter coming from SIFT is used to extract orientation angle .Firstly, the gradient magnitude ( , ) and orientation ( , ) are the bases for the calculation of orientation .They can be calculated utilizing the features from the keypoint and its neighborhoods, as follows [23]: The fundamental principle of the feature tracking method is that persistent surface features common to sequential images can be tracked and identified.The corresponding spatial displacements of these features indicate the ice velocity.The procedures of the rotation-invariant feature-tracking method are described in detail as follows: 1.
Preprocessing.The master image and slave image are strictly georeferenced and orthorectified.By performing this step, the relief displacement and geometric distortion can be removed.

2.
Keypoint localization.In master image, keypoints were chosen at regular intervals.Since the resolution of the MEaSUREs InSAR-Based Antarctica Ice Velocity Map used in this research is 450 m, the corresponding resolution of the Landsat series images is 30 m.Therefore, the interval is 13 pixels.

3.
Orientation assignment.The rotation-invariant parameter coming from SIFT is used to extract orientation angle θ.Firstly, the gradient magnitude m(i, j) and orientation θ(i, j) are the bases for the calculation of orientation θ.They can be calculated utilizing the features from the keypoint and its neighborhoods, as follows [23]: where f (i, j) is the feature for each pixel.In this research, the spectrum is used.Therefore, each f (i, j) is a six-dimensional vector.(i, j) is the locator for each pixel related to the keypoint in the correlation window.The orientation histogram can then be built from these gradient magnitudes and orientations.This orientation histogram has 36 bins covering the 360 degree range of orientations.Each bin is the sum of the gradient magnitudes from the same orientation.The gradient magnitudes for each pixel presented in the histogram are weighted by a Gaussian-weighted circular window.The diameter of the window is five pixels [23], which is an empirical value.Finally, a gradient orientation with respect to the highest bin in the histogram is assigned to the keypoints.4.
Feature tracking.Firstly, reference chip M is rotated to M with orientation θ to eliminate orientation-based displacement (Figure 5).Then, for a candidate feature of the keypoint in the master image, NCC is undertaken on the reference chip centered at the keypoint.According to the rule in [29], the size of the reference chip is 50 × 50, and the corresponding size of the match chip is 100 × 100.Finally, the similarity between M and S is evaluated using the NCC method.5.
Offset calculation.The spatial displacement d can be calculated from the positional difference of conjugate points.Given the known time interval between two image acquisitions, the ice velocity (including flow speed and direction) can be calculated.
Remote Sens. 2018, 10, 42 6 of 16 where ( , ) is the feature for each pixel.In this research, the spectrum is used.Therefore, each ( , ) is a six-dimensional vector.( , ) is the locator for each pixel related to the keypoint in the correlation window.The orientation histogram can then be built from these gradient magnitudes and orientations.This orientation histogram has 36 bins covering the 360 degree range of orientations.Each bin is the sum of the gradient magnitudes from the same orientation.
The gradient magnitudes for each pixel presented in the histogram are weighted by a Gaussianweighted circular window.The diameter of the window is five pixels [23], which is an empirical value.Finally, a gradient orientation with respect to the highest bin in the histogram is assigned to the keypoints.
4. Feature tracking.Firstly, reference chip M is rotated to Mʹ with orientation to eliminate orientation-based displacement (Figure 5).Then, for a candidate feature of the keypoint in the master image, NCC is undertaken on the reference chip centered at the keypoint.According to the rule in [29], the size of the reference chip is 50 × 50, and the corresponding size of the match chip is 100 × 100.Finally, the similarity between Mʹ and S is evaluated using the NCC method.5. Offset calculation.The spatial displacement ∆ can be calculated from the positional difference of conjugate points.Given the known time interval between two image acquisitions, the ice velocity (including flow speed and direction) can be calculated.

Simulated-Image-Based Validation
A simulated image was used to validate the effectiveness of the rotation-invariant featuretracking method.The reason for using a simulated image is due to a lack of in-situ ice velocity data.The simulated image (Figure 6a) as the master image was a 300 × 300 pixel subscene of a Landsat-7 ETM+ image.This image was reprojected into a larger zero-filled grid for image rotation.Figure 6b presents the slave image, which was generated based on the master image by contrarotating by 15°. Figure 6c,d presents the true vector displacement and the result obtained using the rotation-invariant feature-tracking method.The reference chip M for each keypoint in this experiment was a 15-pixel diameter circular window.
The vector displacement maps in Figure 6c,d are smaller than the simulated images, due to the 15-pixel interval from the correlation window of the keypoint.The total number of matching points in the experiment was 256.According to the result shown in Figure 6d, the vector displacement was partially derived using the rotation-invariant feature-tracking method.To further analyze the accuracy of the rotation-invariant feature-tracking method, the number of matching points for different degrees (0°-30°) of rotation is presented in Figure 7. Generally speaking, the rotation degree of a curved outlet glacier will be far less than 30°, so 30° was set as the upper limit of rotation.

Simulated-Image-Based Validation
A simulated image was used to validate the effectiveness of the rotation-invariant feature-tracking method.The reason for using a simulated image is due to a lack of in-situ ice velocity data.The simulated image (Figure 6a) as the master image was a 300 × 300 pixel subscene of a Landsat-7 ETM+ image.This image was reprojected into a larger zero-filled grid for image rotation.Figure 6b presents the slave image, which was generated based on the master image by contrarotating by 15 • .Figure 6c,d presents the true vector displacement and the result obtained using the rotation-invariant feature-tracking method.The reference chip M for each keypoint in this experiment was a 15-pixel diameter circular window.
The vector displacement maps in Figure 6c,d are smaller than the simulated images, due to the 15-pixel interval from the correlation window of the keypoint.The total number of matching points in the experiment was 256.According to the result shown in Figure 6d, the vector displacement was partially derived using the rotation-invariant feature-tracking method.To further analyze the accuracy of the rotation-invariant feature-tracking method, the number of matching points for different degrees (0 • -30 • ) of rotation is presented in Figure 7. Generally speaking, the rotation degree of a curved outlet glacier will be far less than 30 • , so 30 • was set as the upper limit of rotation.When the degree of rotation is less than 5° (especially when it is less than 3°) in Figure 7a, the traditional NCC-based feature-tracking method returns more matches.When the degree of rotation is more than 5°, the matches returned by the traditional NCC-based feature-tracking method rapidly decrease.When the degree of rotation is more than 13°, no match is returned.This indicates that the traditional method is sensitive to the degree of rotation.In contrast, the proposed method is stable at all the different degrees of rotation (Figure 7b).However, when the degree of rotation is less than 5°, the proposed method performs poorly.This is mostly due to the resampling after rotation according to the assigned orientation.Therefore, in this research, the ice velocities for the years are firstly derived by the traditional NCC-based feature-tracking method.And then, the areas in which no velocities are generated are measured using the rotation-invariant feature-tracking method.

Landsat-Image-Based Validation
Only ice velocities of the PRG in 2007, 2008, and 2009 were derived by the different methods from Landsat images, since the majority of the MEaSUREs InSAR-Based Antarctic Ice Velocity Map was acquired from 2007 to 2009. Figure 8 presents the results.
Figure 8a presents profiles from A to B and from C to D, which were chosen to compare the ice velocities.Figure 8b presents the multi-temporal ice velocities from A to B measured by the traditional NCC-based feature-tracking method, the rotation-invariant feature-tracking method, and the velocities from the MEaSUREs InSAR-based Antarctic ice velocity map.The corresponding accuracies are presented in Table 2. Figure 8c presents the multi-temporal ice velocities from C to D measured by the traditional and the improved methods, and the velocities from the MEaSUREs InSAR-based Antarctic ice velocity map.The corresponding accuracies are presented in Table 3.When the degree of rotation is less than 5° (especially when it is less than 3°) in Figure 7a, the traditional NCC-based feature-tracking method returns more matches.When the degree of rotation is more than 5°, the matches returned by the traditional NCC-based feature-tracking method rapidly decrease.When the degree of rotation is more than 13°, no match is returned.This indicates that the traditional method is sensitive to the degree of rotation.In contrast, the proposed method is stable at all the different degrees of rotation (Figure 7b).However, when the degree of rotation is less than 5°, the proposed method performs poorly.This is mostly due to the resampling after rotation according to the assigned orientation.Therefore, in this research, the ice velocities for the years are firstly derived by the traditional NCC-based feature-tracking method.And then, the areas in which no velocities are generated are measured using the rotation-invariant feature-tracking method.

Landsat-Image-Based Validation
Only ice velocities of the PRG in 2007, 2008, and 2009 were derived by the different methods from Landsat images, since the majority of the MEaSUREs InSAR-Based Antarctic Ice Velocity Map was acquired from 2007 to 2009. Figure 8 presents the results.
Figure 8a presents profiles from A to B and from C to D, which were chosen to compare the ice velocities.Figure 8b presents the multi-temporal ice velocities from A to B measured by the traditional NCC-based feature-tracking method, the rotation-invariant feature-tracking method, and the velocities from the MEaSUREs InSAR-based Antarctic ice velocity map.The corresponding accuracies are presented in Table 2. Figure 8c presents the multi-temporal ice velocities from C to D measured by the traditional and the improved methods, and the velocities from the MEaSUREs When the degree of rotation is less than 5 • (especially when it is less than 3 • ) in Figure 7a, the traditional NCC-based feature-tracking method returns more matches.When the degree of rotation is more than 5 • , the matches returned by the traditional NCC-based feature-tracking method rapidly decrease.When the degree of rotation is more than 13 • , no match is returned.This indicates that the traditional method is sensitive to the degree of rotation.In contrast, the proposed method is stable at all the different degrees of rotation (Figure 7b).However, when the degree of rotation is less than 5 • , the proposed method performs poorly.This is mostly due to the resampling after rotation according to the assigned orientation.Therefore, in this research, the ice velocities for the years are firstly derived by the traditional NCC-based feature-tracking method.And then, the areas in which no velocities are generated are measured using the rotation-invariant feature-tracking method.

Landsat-Image-Based Validation
Only ice velocities of the PRG in 2007, 2008, and 2009 were derived by the different methods from Landsat images, since the majority of the MEaSUREs InSAR-Based Antarctic Ice Velocity Map was acquired from 2007 to 2009. Figure 8 presents the results.Figure 8a presents profiles from A to B and from C to D, which were chosen to compare the ice velocities.Figure 8b presents the multi-temporal ice velocities from A to B measured by the traditional NCC-based feature-tracking method, the rotation-invariant feature-tracking method, and the velocities from the MEaSUREs InSAR-based Antarctic ice velocity map.The corresponding accuracies are presented in Table 2. Figure 8c presents the multi-temporal ice velocities from C to D measured by the traditional and the improved methods, and the velocities from the MEaSUREs InSAR-based Antarctic ice velocity map.The corresponding accuracies are presented in Table 3.In Table 3, the rotation-invariant feature-tracking method performed better in 2008 and 2009.In 2007, both the NCC-based feature-tracking method and the rotation-invariant feature-tracking method performed well.Therefore, combining Figure 2 and Tables 2 and 3, it can be concluded that the rotation-invariant feature-tracking method can extend the mapping capability to more complex areas of glacier flow, and can also achieve a high accuracy.Since the computational complexity of the rotation-invariant feature-tracking method is several times higher than the traditional method, the two methods can be combined to derive the ice velocity.Firstly, the traditional NCC-based feature-tracking method is used.Then, the rotation-invariant feature-tracking method measures the ice velocities which are not generated using the traditional method.An example is shown in Figure 2d.

Ice Velocity Variability
A long-term study of the ice velocity of the PRG was undertaken by applying the rotation-invariant feature-tracking method.The ice velocity maps for each year (from 2004 to 2015) are shown in Figure 9.
In Figure 9, the base map for each result comes from a Landsat image in 2004.Due to the Landsat image incompleteness and iceberg calving (between 2013 and 2014), the maps in Figure 9i,j can only present partial ice velocities.Since the accuracy of the ice velocity in 2006 is particularly low, which is possibly due to the snow cover, the ice velocities from 2006 are excluded in the sequential analysis.The changes of the iceberg in front of the PRG are shown in Figure 10.For further analysis, the annual frontal margins and ice velocities are shown in Figures 11 and 12.
The multi-temporal velocities of the iceberg calved between 2013 and 2014 are presented in Figure 12.Point C became open water after 2009 due to the iceberg motion.
To further analyze the multi-temporal ice velocity changes, the elevation change was also taken into consideration.Elevation is a potential influencing factor for ice velocity.An increase in elevation gradient will lead to an increase in ice velocity, while a decrease in elevation gradient will lead to a decrease in ice velocity [30].In this research, the new, enhanced, repeat-track method developed by Hwang et al. [31] was used to calculate the multi-temporal elevations of the upstream part of the PRG from Envisat Sensor Geophysical Data Record (SGDR) data and CryoSat-2 data.The resolution for the bin spacing was about 0.08333 • .The monthly elevation time series was obtained from June 2002 to October 2010 for Envisat and from June 2010 to May 2016 for CryoSat-2.The annual elevations were then obtained by averaging the monthly results.The annual elevations are presented in Figure 13.In Figure 9, the base map for each result comes from a Landsat image in 2004.Due to the Landsat image incompleteness and iceberg calving (between 2013 and 2014), the maps in Figure 9i,j can only present partial ice velocities.Since the accuracy of the ice velocity in 2006 is particularly low, which is possibly due to the snow cover, the ice velocities from 2006 are excluded in the sequential analysis.The changes of the iceberg in front of the PRG are shown in Figure 10.For further analysis, the annual frontal margins and ice velocities are shown in Figures 11 and 12    To further analyze the multi-temporal ice velocity changes, the elevation change was also taken into consideration.Elevation is a potential influencing factor for ice velocity.An increase in elevation gradient will lead to an increase in ice velocity, while a decrease in elevation gradient will lead to a decrease in ice velocity [30].In this research, the new, enhanced, repeat-track method developed by Hwang et al. [31] was used to calculate the multi-temporal elevations of the upstream part of the PRG from Envisat Sensor Geophysical Data Record (SGDR) data and CryoSat-2 data.The resolution for the bin spacing was about 0.08333°.The monthly elevation time series was obtained from June 2002 to October 2010 for Envisat and from June 2010 to May 2016 for CryoSat-2.The annual elevations were then obtained by averaging the monthly results.The annual elevations are presented in Figure 13.

Discussion
We can conclude from Figure 9 that the ice velocities in the PRG vary each year.The glacier and iceberg simultaneously present different ice velocities in their eastern and western areas.Figure 10 indicates that the iceberg in front of the PRG calved into pieces between 2013 and 2014, which then drifted away.Meanwhile, the iceberg in front of the Polar Time Glacier also calved during this period; however, it remained in the Prydz Bay area after calving.High velocities in the frontal area of the Polar Times Glacier in 2008 and 2012 result from big movement of its frontal iceberg, which may lead To further analyze the multi-temporal ice velocity changes, the elevation change was also taken into consideration.Elevation is a potential influencing factor for ice velocity.An increase in elevation gradient will lead to an increase in ice velocity, while a decrease in elevation gradient will lead to a decrease in ice velocity [30].In this research, the new, enhanced, repeat-track method developed by Hwang et al. [31] was used to calculate the multi-temporal elevations of the upstream part of the PRG from Envisat Sensor Geophysical Data Record (SGDR) data and CryoSat-2 data.The resolution for the bin spacing was about 0.08333°.The monthly elevation time series was obtained from June 2002 to October 2010 for Envisat and from June 2010 to May 2016 for CryoSat-2.The annual elevations were then obtained by averaging the monthly results.The annual elevations are presented in Figure 13.

Discussion
We can conclude from Figure 9 that the ice velocities in the PRG vary each year.The glacier and iceberg simultaneously present different ice velocities in their eastern and western areas.

Discussion
We can conclude from Figure 9 that the ice velocities in the PRG vary each year.The glacier and iceberg simultaneously present different ice velocities in their eastern and western areas.Figure 10 indicates that the iceberg in front of the PRG calved into pieces between 2013 and 2014, which then drifted away.Meanwhile, the iceberg in front of the Polar Time Glacier also calved during this period; however, it remained in the Prydz Bay area after calving.High velocities in the frontal area of the Polar Times Glacier in 2008 and 2012 result from big movement of its frontal iceberg, which may lead by the interaction between the Polar Times Glacier and the Polarforschung Glacier.The ice velocity in 2008 was generated using a pair of Landsat images (December 2007 and March 2009) whose time interval is 432 days including two summers.In general, the glacier presents the highest flow speed in summer.Therefore, this may be the other reason leading to the high flow speed in 2008.The interaction between the glacier terminus and iceberg can hamper the movement of the glacier tongue.Figure 14 presents the gaps between the glacier terminus and iceberg in 2005 and 2006.It can be found that the outlet glacier moved faster than the iceberg during this period.Sea ice changes can influence the motion of the outlet glacier and iceberg [11].The lower speed of the frontal iceberg of the PRG may result from the surrounded sea ice.Therefore, more multi-source satellite images are needed to prove this.
Remote Sens. 2018, 10, 42 13 of 16 in summer.Therefore, this may be the other reason leading to the high flow speed in 2008.The interaction between the glacier terminus and iceberg can hamper the movement of the glacier tongue.Figure 14 presents the gaps between the glacier terminus and iceberg in 2005 and 2006.It can be found that the outlet glacier moved faster than the iceberg during this period.Sea ice changes can influence the motion of the outlet glacier and iceberg [11].The lower speed of the frontal iceberg of the PRG may result from the surrounded sea ice.Therefore, more multi-source satellite images are needed to prove this.
(a) (b) Figure 11 presents annual frontal margins and multi-temporal ice velocities of the PRG through nine uniformly distributed points.In Figure 11b, it is clear that the ice flow direction of the PRG is not straight, but to the northeast, and the ice flow direction of the Polar Times Glacier located on the PRG's west flank is to the northwest.This can likely be attributed to the interaction of the Polarforschung Glacier and the Polar Times Glacier. Figure 11c shows that the variation trends of the ice velocities in the PRG are similar at each point, but the interannual change varies.The ice velocities in 2010 and 2011 were relatively low.In contrast, the ice velocities in 2008 and 2012 were higher.Linear regression was used to analyze this annual change for each point.The R 2 values are all less than 0.05.This indicates that there was no significant linear trend for the ice velocities of the PRG during this period.
Several particular characteristics can also be observed in Figure 11c.The nine points can be separated into three groups according to their longitudes, i.e., ADG, BEH, and CFI.The lowest ice velocities are observed in ADG; the highest ice velocities are observed in CFI; and the ice velocities of BEH are between these two groups.This is possibly a result of the PRG rotation.The flow direction of the PRG is to the northeast, which results in a higher velocity in the west side of the PRG and a lower velocity in the east side.In addition, the nine points can also be separated into three groups according to their latitudes, i.e., ABC, DEF, and GHI.The ice velocities of ABC are the lowest, and the ice velocities of GHI are the highest.The ice velocity of each point gradually increases from the south to the north.However, the ice velocity of point F does not abide by this rule.This may suggest that the stress coming from the interaction between the PRG and the Polar Times Glacier plays a more important role in the ice velocity at point F than the stress coming from the upper glacier.
Although the PRG calved between 1973 and 1989, the calved iceberg did not move away from the terminus.This results in a similar variation trend in their ice velocities from 2004 to 2012, which is clearly presented in Figures 11c and 12d.However, due to the annual change of the contact points between the PRG and the calved iceberg (Figure 12c), the change trends are slightly different, especially in 2011 and 2012.This suggests that the motion of the iceberg mainly originates from the PRG.Combining Figure 12b,d, the velocities of each point are highly consistent between 2004 and 2008.This suggests that the contact points and motion directions did not change during this period.After 2008, the changes in the velocity of each point are possibly due to the changing of the contact points, as mentioned before.Figure 11 presents annual frontal margins and multi-temporal ice velocities of the PRG through nine uniformly distributed points.In Figure 11b, it is clear that the ice flow direction of the PRG is not straight, but to the northeast, and the ice flow direction of the Polar Times Glacier located on the PRG's west flank is to the northwest.This can likely be attributed to the interaction of the Polarforschung Glacier and the Polar Times Glacier. Figure 11c shows that the variation trends of the ice velocities in the PRG are similar at each point, but the interannual change varies.The ice velocities in 2010 and 2011 were relatively low.In contrast, the ice velocities in 2008 and 2012 were higher.Linear regression was used to analyze this annual change for each point.The R 2 values are all less than 0.05.This indicates that there was no significant linear trend for the ice velocities of the PRG during this period.
Several particular characteristics can also be observed in Figure 11c.The nine points can be separated into three groups according to their longitudes, i.e., ADG, BEH, and CFI.The lowest ice velocities are observed in ADG; the highest ice velocities are observed in CFI; and the ice velocities of BEH are between these two groups.This is possibly a result of the PRG rotation.The flow direction of the PRG is to the northeast, which results in a higher velocity in the west side of the PRG and a lower velocity in the east side.In addition, the nine points can also be separated into three groups according to their latitudes, i.e., ABC, DEF, and GHI.The ice velocities of ABC are the lowest, and the ice velocities of GHI are the highest.The ice velocity of each point gradually increases from the south to the north.However, the ice velocity of point F does not abide by this rule.This may suggest that the stress coming from the interaction between the PRG and the Polar Times Glacier plays a more important role in the ice velocity at point F than the stress coming from the upper glacier.
Although the PRG calved between 1973 and 1989, the calved iceberg did not move away from the terminus.This results in a similar variation trend in their ice velocities from 2004 to 2012, which is clearly presented in Figures 11c and 12d.However, due to the annual change of the contact points between the PRG and the calved iceberg (Figure 12c), the change trends are slightly different, especially in 2011 and 2012.This suggests that the motion of the iceberg mainly originates from the PRG.Combining Figure 12b,d, the velocities of each point are highly consistent between 2004 and 2008.This suggests that the contact points and motion directions did not change during this period.After 2008, the changes in the velocity of each point are possibly due to the changing of the contact points, as mentioned before.
Combining Figures 11 and 13, it can be clearly seen that there was no significant change in elevation and ice velocity for the PRG from 2004 to 2015.The elevation decreased by ~2.3 m, corresponding to an ice velocity increase of ~10 m (mean value for each point in Figure 11a).Except for two periods (from 2009 to 2010 and from 2014 to 2015), the interannual changes of elevation and ice velocity were also consistent.With the elevation increasing in three periods (from 2007 to 2008, from 2011 to 2012, from 2013 to 2014), the ice velocities also increased.With the elevation decreasing in two periods (from 2008 to 2009, from 2012 to 2013), the ice velocities decreased.This phenomenon proves the statement that elevation is a potential driving factor for ice velocity, as mentioned in [28].However, an opposite trend is shown between the elevation and ice velocity of the PRG from 2009 to 2010.By combining Figure 9f, a possible reason for this may be that the obvious deflection of the PRG to the northeast led to low velocities in the middle and east side of the PRG, which made the forward velocity decrease.From 2014 to 2015, the opposite trend between the elevation and ice velocity of the PRG can be attributed to the calving of the iceberg.

Conclusions
Ice velocity is a fundamental characteristic of glacier dynamics.In this research, annual ice velocity variations for the PRG between 2004 and 2015 were investigated from sequential Landsat images using a feature tracking method.In order to overcome the shortcoming of the traditional NCC-based feature-tracking method, which cannot measure ice velocities in some areas where the glacier curves, a rotation-invariant parameter was introduced, and a rotation-invariant feature-tracking method was creatively built.
The validations were implemented based on a simulated image and a Landsat image.The results with the simulated image reveal that the rotation-invariant feature-tracking method can measure ice velocities from more areas of the outlet glacier, especially when the rotation angle of the glacier is more than 5 • .Moreover, the proposed method performs well in the ice velocity comparison between the traditional NCC-based feature-tracking method, the rotation-invariant feature-tracking method, and the velocities from the MEaSUREs InSAR-Based Antarctic Ice Velocity Map.
Ice velocities in the PRG have varied each year.However, the trend from 2004 to 2015 shows no significant change, in that the ice velocity increased by only about 10 m during this period.The maximum ice velocity of the frontal margin of the PRG reached about 900 m/a, and the maximum ice velocity of the iceberg reached about 1000 m/a before calving.The interaction of the Polar Times Glacier and the Polarforschung Glacier impacted the direction and the displacement of the PRG.This impact also led to a higher velocity in the west side of the PRG and a lower velocity in the east side.Elevation changes partly influenced the ice velocities of the PRG; however, this impact was easily overshadowed by the iceberg calving.

Figure 2 .
Figure 2. Ice velocity maps for the Polar Times Glacier.(a) Master image from February 2005; (b) Slave image from February 2006; (c) Ice velocity results obtained using the traditional NCC-based featuretracking method; (d) Ice velocity results obtained combing the traditional NCC-based feature-tracking method, the rotation-invariant feature-tracking method.

Figure 2 .
Figure 2. Ice velocity maps for the Polar Times Glacier.(a) Master image from February 2005; (b) Slave image from February 2006; (c) Ice velocity results obtained using the traditional NCC-based feature-tracking method; (d) Ice velocity results obtained combing the traditional NCC-based feature-tracking method, the rotation-invariant feature-tracking method.

Figure 3 .
Figure 3. Schematic diagram for the pixel-level ice flow of a curved outlet glacier.M is the reference chip.S and Sʹ are match chips.The red pixel in the slave image is the matching point which corresponds to the keypoint in the master image.

Figure 4 .
Figure 4.The flowchart of the rotation-invariant feature-tracking method.

Figure 3 .
Figure 3. Schematic diagram for the pixel-level ice flow of a curved outlet glacier.M is the reference chip.S and S are match chips.The red pixel in the slave image is the matching point which corresponds to the keypoint in the master image.

Figure 3 .
Figure 3. Schematic diagram for the pixel-level ice flow of a curved outlet glacier.M is the reference chip.S and Sʹ are match chips.The red pixel in the slave image is the matching point which corresponds to the keypoint in the master image.

Figure 4 .
Figure 4.The flowchart of the rotation-invariant feature-tracking method.

Figure 4 .
Figure 4.The flowchart of the rotation-invariant feature-tracking method.

Figure 5 .
Figure 5. Schematic diagram for the reference chip rotation.M is the reference chip.Mʹ is the rotated chip.S is the match chip.The red pixel in the slave image is the matching point which corresponds to the keypoint in the master image.

Figure 5 .
Figure 5. Schematic diagram for the reference chip rotation.M is the reference chip.M is the rotated chip.S is the match chip.The red pixel in the slave image is the matching point which corresponds to the keypoint in the master image.

Figure 6 .Figure 7 .
Figure 6.Ice velocity validation based on the simulated image.(a) Simulated master image; (b) Simulated slave image after 15° manual contrarotation based on the master image; (c) True vector displacement map; (d) Vector displacement derived using the rotation-invariant feature-tracking method.

Figure 6 .Figure 6 .Figure 7 .
Figure 6.Ice velocity validation based on the simulated image.(a) Simulated master image; (b) Simulated slave image after 15 • manual contrarotation based on the master image; (c) True vector displacement map; (d) Vector displacement derived using the rotation-invariant feature-tracking method.

Figure 7 .
Figure 7. Number of matching points for different degrees of rotation (0-30 • ).(a) Result of the traditional NCC-based feature-tracking method; (b) Result of the rotation-invariant feature-tracking method.

Figure 8 .
Figure 8. Ice velocity comparison between the traditional NCC-based feature-tracking method, the rotation-invariant feature-tracking method, and the velocities from the MEaSUREs InSAR-based Antarctic ice velocity map.(a) The location of the profiles from A to B and from C to D; (b) The different ice velocities from A to B; (c) The different ice velocities from C to D.

Figure 8 .
Figure 8. Ice velocity comparison between the traditional NCC-based feature-tracking method, the rotation-invariant feature-tracking method, and the velocities from the MEaSUREs InSAR-based Antarctic ice velocity map.(a) The location of the profiles from A to B and from C to D; (b) The different ice velocities from A to B; (c) The different ice velocities from C to D.

Figure 10 .
Figure 10.Sequential Landsat series images presenting changes of the iceberg in front of the PRG.The glacier outlines shown in different colors except November 2013.Blue line area presents the outline of the frontal iceberg in January 2013.(a) January 2013; (b) November 2013; (c) November 2014; (d) April 2016.

Figure 10 .Figure 11 .
Figure 10.Sequential Landsat series images presenting changes of the iceberg in front of the PRG.The glacier outlines shown in different colors except November 2013.Blue line area presents the outline of the frontal iceberg in January 2013.(a) January 2013; (b) November 2013; (c) November 2014; (d) April 2016.Remote Sens. 2018, 10, 42 11 of 16

Figure 11 .Figure 12 .
Figure 11.Annual frontal margins and multi-temporal ice velocities of the PRG.(a) Nine uniformly distributed points named "A" to "I" were chosen on the PRG for further analysis; (b) Annual frontal margins from 2004 to 2015 shown in different colors; (c) Annual ice velocities of different points from 2004 to 2015.

Figure 13 .
Figure 13.Annual elevation results from Envisat and CryoSat-2 (2002-2016).The blue line shows the elevation from Envisat (June 2002 to October 2010), and the orange line shows the elevation from CryoSat-2 (June 2010 to May 2016).

Figure 12 .Figure 12 .
Figure 12.Annual contour lines and multi-temporal velocities of the calved iceberg.(a) Nine uniformly distributed points named "A" to "I" were chosen on the PRG; (b) Annual contour lines shown in different colors; (c) Annual contact points of the PRG and the calved iceberg shown in different colors; (d) Annual ice velocities of the different points.

Figure 13 .
Figure 13.Annual elevation results from Envisat and CryoSat-2 (2002-2016).The blue line shows the elevation from Envisat (June 2002 to October 2010), and the orange line shows the elevation from CryoSat-2 (June 2010 to May 2016).
Figure 10 indicates that the iceberg in front of the PRG calved into pieces between 2013 and 2014, which then drifted away.Meanwhile, the iceberg in front of the Polar Time Glacier also calved during this period; however, it remained in the Prydz Bay area after calving.High velocities in the frontal area of the Polar Times Glacier in 2008 and 2012 result from big movement of its frontal iceberg, which may lead by the interaction between the Polar Times Glacier and the Polarforschung Glacier.The ice velocity in 2008 was generated using a pair of Landsat images (December 2007 and March 2009) whose time interval is 432 days including two summers.In general, the glacier presents the highest flow speed

Figure 13 .
Figure 13.Annual elevation results from Envisat and CryoSat-2 (2002-2016).The blue line shows the elevation from Envisat (June 2002 to October 2010), and the orange line shows the elevation from CryoSat-2 (June 2010 to May 2016).

Table 1 .
Sequential Landsat series images of the PRG.

Table 1 .
Sequential Landsat series images of the PRG.

Table 2 .
The accuracies of the different methods compared with the MEaSUREs InSAR-based Antarctic ice velocity map (m/a) (Profile from A to B).

Table 3 .
The accuracies of the different methods compared with the MEaSUREs InSAR-based Antarctic ice velocity map (m/a) (Profile from C to D).

Table 2
, the RMSEs for 2009 are far greater than the RMSEs in 2007 and 2008 compared with the MEaSUREs InSAR-based Antarctic ice velocity map.The lower the RMSE, the better the satellitebased ice velocity.The ice velocities in 2007 and 2008 are more accurate than those in 2009.The peak

Table 2 .
The accuracies of the different methods compared with the MEaSUREs InSAR-based Antarctic ice velocity map (m/a) (Profile from A to B).

Table 3 .
The accuracies of the different methods compared with the MEaSUREs InSAR-based Antarctic ice velocity map (m/a) (Profile from C to D).

Table 2 ,
the RMSEs for 2009 are far greater than the RMSEs in 2007 and 2008 compared with the MEaSUREs InSAR-based Antarctic ice velocity map.The lower the RMSE, the better the satellite-based ice velocity.The ice velocities in 2007 and 2008 are more accurate than those in 2009.The peak signal-to-noise ratio (PSNR) is also used in this validation.The higher the PSNR, the better the satellite-based ice velocity.The results of the PSNRs also indicate that the ice velocities in 2007 and 2008 are more accurate.Comparing RMSE and PSNR, it can be found that the NCC-based feature-tracking method performed better in 2007, but worse in 2008 and 2009.In contrast, the rotation-invariant feature-tracking method performed better in 2008 and 2009, but worse in 2007.