Improved multi-scale image matching approach for monitoring Amery ice shelf velocity using Landsat 8

ABSTRACT The velocity of an ice shelf is an important characteristic to understand its dynamics and interaction with the internal ice sheet. Therefore, we develop an improved multi-scale image matching method for producing a complete and accurate Amery ice shelf velocity field from Landsat 8 images. First, we investigate the relationship between the template size and the image entropy and propose an image entropy-based preliminary operation for distinguishing the high-contrast regions from the low-contrast regions prior to iteratively determining the optimum template size. Second, a Gaussian pyramid image-based hierarchical data structure is designed to support a coarse-to-fine image matching strategy.The image entropy-based matching method is applied on the top layer of the image pyramid to guarantee the matching results have optimal completeness and high accuracy. Finally, a postprocess procedure is performed to derive a complete and accurate Amery ice shelf velocity field. Experimental results demonstrate that the proposed method can significantly improve computational efficiency. Moreover, the proposed method provides more accurate and robust matching results than other existing methods, particularly over the low-contrast surfaces. Additionally, the derived velocity field also shows good consistency with the one acquired from MEaSUREs Annual Antarctic Ice Velocity Maps 2016–2017.


Introduction
Antarctic ice shelves are significant components of ice sheets due to the ice-ocean-atmosphere interface and vulnerability to atmospheric and oceanic temperature changes at the regional and global scale (Mercer, 1978;Vaughan & Doake, 1996). Compared with groundbased ice sheets or continental glaciers, ice shelves are more sensitive to climate changes (Chen et al., 2013). Mapping the ice velocities of ice shelves remarkably enhanced our understanding of the physics of ice flow, the mass balance of ice sheets and their contribution to sea level (Mouginot & Rignot, 2017).
To date, many automatic approaches for monitoring the ice shelf velocity from remotely sensed data have been developed, such as Synthetic Aperture Radar (SAR) image-based interferometry, SAR image-based speckle tracking and satellite visible-infrared imagebased feature tracking. Among these existing methods, the application of SAR interferometry achieved the highest accuracy. However, it is not flexible to unwrap the unrecoverable phases, especially in fast-moving areas (Lee, Seo, Han, Yu, & Scambos, 2012). In contrast, the SAR image-based speckle tracking determines the displacements by using the cross-correlation method of the speckle patterns (Joughin, 2002), which can effectively address the limitation of interferometry. As a result, the combined methods of interferometry and speckle tracking generated the first ice velocity maps at the ice sheet scale (Jezek, Sohn, & Noltimier, 1998;Joughin, Smith, Howat, Scambos, & Moon, 2010;Rignot, Mouginot, & Scheuchl, 2011b). However, there were limited available free SAR data sources. The satellite visible-infrared image-based feature tracking, especially the use of Landsat 8 satellite images is increasingly popular due to its higher radiometric quantization and approximately half-pixel geolocation accuracy which made it possible to track subtle surface features (Jeong & Howat, 2015), such as sastrugi.
The surface velocity of ice shelves can be formulated into the displacement of corresponding pixels in a single image pair over a time span. To obtain these displacements, various image matching-based methods were tested and applied for feature tracking beyond repeated images. Berthier, Raup and Scambos (2003) and Kääb (2002) employed normalized cross-correlation (NCC) algorithm in spatial domain to determine a velocity field. Kaufmann and Ladstädter (2003) adopted least squares matching for measuring surface deformation and flow velocity of rock glaciers using digital photogrammetry. Compared to those in the spatial domain, the image matching methods in the frequency domain, as a convolution operation, is less time-consuming (Heid & Kääb, 2012). Cecilie, Jostein, Jon-Ove and Bengt (1997) calculated the cross-correlation in frequency domain by utilizing fast Fourier transform. Some researchers (Leprince, Barbot, Ayoub, & Avouac, 2007;Quincey & Glasser, 2009;Scherler, Leprince, & Strecker, 2008) also applied the phase correlation approach to obtain ice velocity. Fitch et al.(2002, September) carried out the orientation correlation algorithm on orientation images for a more robust ice velocity field.
Although many studies have been done in recent years, automatically mapping velocity field of ice shelf from remotely sensed data accurately, completely and efficiently remains a great challenge because of the following reasons. First, repeated image matching techniques are widely used to track persistent features on the ice surface for deriving the ice motion information. As Ahn and Howat (2011) have noted, the performance of the conventional image matching method depends on the template size. Debella-Gilo and  calculated the locally optimum template size by finding the first peak of NCC coefficients from a series of template sizes. However, the way to find the locally adaptive template always results in numerous computations. Second, the identification of the matchable templates and their presences (or absences) in the search image is of significance for reducing the occurrences of the bad and spurious match. As a matter of fact, a good matching candidate over an area is easily obtained if the surface has good contrast (such as crevasses) and the images contain low level of noise. However, highly reflective targets, such as snow and sastrugi usually contain limited radiometric information over ice, which poses a challenge for obtaining an exact match with precise location.
To address these above issues, this paper develops an improved multi-scale image matching method. Firstly, we investigate the relationship between the template size and the image entropy and propose an image entropy-based preliminary operation for distinguishing the high-contrast regions from the low-contrast regions before determining the optimum template size. Moreover, to reduce the computational complexity and enhance the matching efficiency, we design a Gaussian pyramid image-based hierarchical data structure to support a coarse-to-fine image matching strategy. Starting from the coarsest image at the top pyramid level, we conduct the improved multi-scale image matching approach and the matching results on the coarse-resolution image as the seeds are progressively propagated and expanded to the higher resolution image to obtain more accurate motion information on the ice surface.
The rest of this paper is organized as follows. The Methodology section describes the proposed method in detail. The Experimentation and analysis section presents the experimental results and analysis for evaluating the proposed method. This paper concludes with a discussion of future research considerations in the Conclutions section.

Methodology
In this paper, we develop an improved multi-scale image matching method for producing a complete and accurate Amery ice shelf velocity field from Landsat 8 images. As illustration in Figure 1, the proposed method consists of the following steps: (1) Preprocessing of original imagery; (2) Improved multi-scale image matching method for generating displacements of the single image pair; (3) Ice shelf velocity field production based on image mosaic and median filtering.

Preprocessing of original imagery
As shown in Figure 1, we need to track surface features between image pairs of Landsat 8 Operational Land Imager (OLI) panchromatic channel (Band 8) with the resolution of 15 Â 15. To improve the quality of the feature tracking algorithm, the standard image enhancement procedure is applied to enhance the surface feature contrasts Glasser & Scambos, 2008), such as crevassing, sastrugi and smaller-amplitude roughness. Histogram equalization is the most commonly used method for its simplicity and relatively good performance on almost all types of images (Abdullah-Al-Wadud, Kabir, Dewan, & Chae, 2007), which performs its operation by calculating the probability distribution of the input level and remapping the output gray level. Figure 2 exhibits the comparison of results between enhanced image and non-enhanced image and their associated histograms, and Figure 3 shows the detailed comparisons of some typical surface features. As shown in Figure 3, the preprocessing of original imagery adjusts the image contrast and enhances the effectiveness of embedded information to better help feature identification for image matching. Consequently, the enhanced images will serve as the input of the feature tracking method.
Improved multi-scale image matching for generating displacements of the single image pair The procedure for multi-scale image matching based on Gaussian pyramid of Landsat 8 images is illustrated in Figure 4. We firstly design the hierarchical data structure based on the Gaussian pyramid imagery, and then distinguish between pixels in high-contrast and low-contrast areas by using image entropy on the coarsest image at the top pyramid level (Reference imagery and Search imagery G r3 ). For different type of areas, different range of template sizes are applied to find an optimum template size. The obtained accurate and complete initial positions are used for the next layer matching. The whole image matching is not accomplished until the matching reaches the last layer of the Gaussian pyramid (G r0 ) which has the highest resolution. Key algorithms of the improved multi-scale image matching approach are given in more detail below.

NCC in the frequency domain
As aforementioned, repeated image matching techniques, have become increasingly popular in analyzing and monitoring ice surface movements in recent years. Among these matching techniques, the areabased method with NCC as its similarity measure is widely used because of its simplicity and reliability. Also, its correlation computation can be speeded up by the fast Fourier transformation (Heid & Kääb, 2012).
For the NCC-based image matching method, the persistent features on the ice surface are tracked to estimate a displacement between the earlier image (also called the reference image) and the later image (also called the search image). More specifically, an image subset (or a template) is first selected in the reference image centered on pixel x; y ð Þ. Then, an image subset of same size is extracted from the search image but translated by x 0 ; y 0 ð Þ pixels inside a predefined searching window, which can be defined based on the maximum possible velocity, the image resolution and time span of both the reference image and the search image (Terzopoulos, 1986), and compared with NCC coefficient that can be calculated using Equation (1) as follows: Where T x; y ð Þ represents the pixel value of the center position of the template window, I x 0 ; y 0 ð Þ represents the pixel value in the search image, and n represents the number of pixels in the template window. The T 0 x; y ð Þ denotes the average of gray values of the pixels in the template window subtracted from the T x; y ð Þ. I 0 x 0 ; y 0 ð Þ represents the average of gray values of the pixels in the search window subtracted from the I x 0 ; y 0 ð Þ. This operation is repeated for different values of x; y ð Þ and the position with the maximum correlation coefficient is considered as the true displacement.

Relationship among template sizes, NCC coefficients and displacements
With regards to the NCC-based image matching method, the template size is a vital parameter for the performance of feature tracking. As Debella-Gilo and Kääb (2012) mentioned, adaptive template sizes can be derived from the first local peak of variations of the NCC coefficients on the same central pixel with different template sizes (as shown in Figure 5). This method is not available for low-contrast areas, such as shadows, snow, and other monotonous surfaces that lack the remarkable signal variability (as shown in Figure 5(d)). However, lacking sufficient features (or low-contrast in this paper) is a common issue on ice shelf.
In our implementation, several typical image positions with sufficient and insufficient signal variability are selected to analyse the relationship among template sizes, NCC coefficients and displacements. Using the older image as reference, the maximum NCC coefficients are calculated for different template sizes ranging from 9 Â 9 to 171 Â 171. As shown in Figure 5(a), for each template size, only the first local peak of the NCC coefficient and the position of the match are recorded, and the first peak is calculated according the previous work (Debella-Gilo & , as defined in Equation (2). Where ρ denotes NCC coefficients, τ denotes template size. It is confirmed that for the area with sufficient features (as shown in Figure 5(b)), the displacements in the x-and y-direction are fixed around the calculated adaptive template. In terms of the area with low-contrast (as shown in Figure 5(d)), the displacements calculated by the adaptive template exceed the prior velocity range (Mouginot & Rignot, 2017) and have a huge fluctuation around the adaptive template. Alternatively, if we apply the bigger initial template sizes the surface features on the other areas will be included, and the calculated displacement is valid by using calculated adaptive template (as shown in Figure 5(e)). Consequently, bigger initial template sizes can be used to find the valid displacement of area with low contrast.
Distinguishing between pixels on high-contrast and low-contrast areas using image entropy As described in previous section, the effectiveness and robustness of local adaptive template size depends on whether signal variations are sufficient or not. It is well known that the surface features are associated with the local grayscale variations in the digital images, which is primarily determined by physical and geometric characteristics of the scene. For the digital images, the existence of features indicates that the area has more information and vice versa. Both signal-to-noise ratio (SNR) (Immerkaer, 1996) and image entropy (Shannon, 1948) have been the powerful tools for measuring the information that the image contains. Image entropy is a measure of the amount of image information. The sufficient surface features mean the area contains more information, and its associated image entropy is bigger (Figure 6(a)). Otherwise, the image entropy is smaller (Figure 6(b)). A typical example about the relationship between SNR and template sizes, and the relationship between image entropy and template sizes for the low-contrast areas and the high-contrast areas is presented in Figure 7. Experimental results suggest that the image entropy of the high-contrast areas always exceeds that of low-contrast areas with the template size increasing. On the contrary, Figure 5. Typical areas and their associated relationships with the template sizes, the maximum NCC coefficients and the displacements: (a) the relationship between the template sizes, the maximum NCC coefficients and the displacements of highcontrast area; (b) area with high-contrast(template size 39 Â 39); (c) the relationship between the template sizes, the maximum NCC coefficients, and the displacements of low-contrast area; (d) area with low-contrast(template sizes 39 Â 39); (e) the relationship between the template sizes, the maximum NCC coefficients and the displacements of low-contrast area (starting from bigger initial template size); (f) area with low-contrast (starting from bigger initial template size(109 Â 109).
when the template size is small, the SNR of lowcontrast area even exceeds that of high-contrast area (Figure 7(c)). As a result, the image entropy is adopted as a measurement for distinguishing the low-contrast areas and the high-contrast areas on the optical images.
As shown in Figure 7, we also conclude that the template size 23 and its associated entropy can be used to distinguish low-contrast pixels from high-contrast pixels. In our implementation, we manually mark some high-contrast and low-contrast areas and the image entropies of the corresponding positions are  calculated. As a result, the positions with image entropy less than 8.9 are marked as the low-contrast areas, while the rests are marked as the high-contrast areas. Regions over Landsat 8 image are marked as the low-contrast regions and the high-contrast regions using the image entropy (as illustrated in Figure 8).
Multi-scale image matching based on Gaussian pyramid using the adaptive template size The computational efficiency of iteratively finding the locally adaptive template size for each matching position needs to be improved. Therefore, to reduce the computational complexity and enhance the matching efficiency, we conduct multi-scale image matching based on Gaussian pyramid (Adelson, Anderson, & Bergen, 1984), which is widely used.
As shown in Figure 4, prior to the NCC-based image matching in the frequency domain, the matching positions on the query reference image G r3 are classified into low-contrast areas and high-contrast areas based on image entropy. After that, we follow the idea of Debella-Gilo and  and determine the locally adaptive template size for each matching position by setting different initial template size for two types of pixels: for high-contrast pixels the range of template sizes is 9 Â 9 to 171 Â 171,while for low-contrast pixels, the range of template sizes is 85 Â 85 to 171 Â 171. Then, the NCC-based image matching in the frequency domain is performed on the coarsest image at the top pyramid level. As a matter of fact, the displacement of x-and y-directions should be fixed around the optimum template size. To improve the matching accuracy and eliminate the erroneous matches, we need to check whether the true match is obtained at the optimum template size. As defined in Equation (3), x 1 , x 2 , y 1 , y 2 correspond to the coordinates of the displacement in the x-and y-directions for two successive template sizes during the process of increasing the template size progressively. If the matching position does not change within that range δ, it is considered as a true match. Otherwise, it is considered as a bad or spurious match and therefore removed. ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi The matching result from the coarsest image at the top pyramid level is then used as initial seeds to search for more accurate matches at the next level of higher resolution image. During the multi-scale image matching procedure, the low-resolution displacements from a coarse-resolution level G i need to be propagated to the next higher resolution level G i−1 in the pyramid. For the inter-level propagation, we employ an inverse distance weighted interpolation approach (Liu, Wang, Tang, & Jezek, 2012) to estimate approximate displacements for each matching position at the next level G i-1 based on the matching results at the query level G i . The iterative process repeats till the matching step is performed at the bottom level in the pyramid. As a result, a coarse-tofine matching result is derived from the bottom level image with the original image resolution.

Postprocessing
After processing all the selected image pairs mentioned above, we derive a complete and accurate velocity field of the Amery ice shelf with a successive scheme, which includes outlier removal, integration into a single velocity and gaps filling (Shannon, 1948).
Outlier removal Due to thin clouds or changing snow conditions between the images (Liu et al., 2012), not all features on the ice shelf surface can be tracked persistently with optical image matching, which probably results in some spurious matches. Thus, we in this section use a neighbor filter to remove these spurious matches based on the assumption of a certain level of minimum smoothness in the ice shelf surface velocity field. In our implementation, we compute the differences between velocity values in the x-and y-direction at the query pixel and its 24 neighbors, respectively. If a velocity value in the x-or y-direction is three times more than the standard deviation Δ (Equation (4)), it is considered as an outlier (Ahn & Howat, 2011).
where v x; y ð Þ represents the velocity in the x-and y-direction of the query pixel, and v x i ; y i ð Þ represents the velocity in the x-and y-direction of neighboring pixels.

Integration into a single velocity and gaps filling
Displacements of the single image pair can be generated using the approach above. To obtain the whole ice velocity field of the Amery ice shelf, in this section, we need to integrate eight image pairs into a single velocity and use median filter for interpolating image gaps. Figure 9(a) shows the overlapping images (image pair of World Reference System 2(WRS2) Path 127 and Row 111, image pair of WRS2 Path 129 and Row 110) from different collection time, and each single image pair will generate a velocity field. Due to the overlapping areas among contiguous Landsat 8 images (as shown in Figure 9(b)), there are multiple different velocities at the overlapping area. Therefore, a key postprocessing step during mosaic is the determination of the best single displacement. In our implementation, we follow the idea of the previous work (Fahnestock et al., 2016) and adopt a weighted average-based method to determine the optimal ice velocities since the weighted average-based method emphasized the best-correlated, longest time-interval data. As defined in Equation (5), the weighted average velocity value for each pixel is determined based on image pair time-interval (time/ 365) multiplied by the square root of the correlation coefficient of NCC method. For the non-overlapping area, we preserve the calculated ice velocity in this area. Finally, we conduct gap filling using 5 × 5 median filter because it is powerful for reducing random noises, especially when the noise amplitude probability density has large tails, and periodic patterns (Huang, Yang, & Tang, 1979). As a result, Moderate resolution Imaging Spectroradiometer (MODIS) Mosaic of Antarctica is used as a reference frame (Scambos, Haran, Fahnestock, Painter, & Bohlander, 2007) to generate the final production of ice velocity.
where W denotes the weight of each pixel, t denotes timeinterval of image pair, R denotes the correlation coefficient, V denotes the weighted average velocity value.

Experimentation and analysis
Amery ice shelf as the third largest ice shelf in the Antarctic (as shown in Figure 10) covers an area of 1.4 × 10 6 km 2 and is the largest drainage basin in East Antarctica with approximately 16% of the ice in this sector of the continent draining through it (Budd, 1966), with flow velocities of 350-1700 m/ year (Allison, 2003). In order to further define the scope of the experimental area, we used the (Making Earth System Data Records for Use in Research Environments) MEaSUREs Antarctic Figure 9. A simple example of overlapping images derived from image mosaic.
Grounding Line (Rignot, Mouginot, Morlighem, Seroussi, & Scheuchl, 2014;Rignot, Mouginot, & Scheuchl, 2011a), the point at which these glaciers start to float is the grounding line. The location of the grounding Line is important because mass balance from Antarctica is strongly associated with variations in the ice shelves and their grounding Lines.

Experimental data
The Landsat 8 on-board OLI sensor launched in February 2013 provides imagery with 12-bit radiometric resolution, a 16-fold increase over Enhanced Thematic Mapper Plus (ETM+). Brunt, Fricker, Padman, Scambos, and O'Neel (2010) investigated and found out that radiometric performance such as SNR and its dynamic range of OLI outperforms those of ETM+. Since Landsat 8 images can be used as good candidates than other optical images, the panchromatic images (band 8, 15 m pixel resolution) were used as experimental data from L1GT imagery, where geometrically corrected products can be systematic terrain corrected (Storey, Choate, & Lee, 2014), which are freely available from United States Geological Survey. The time-intervals of the experimental images in this paper are around 1 year, in which the longest interval is 368 days and the shortest interval is 352 days.

Experimental analysis
Effectiveness of distinguishing the low-contrast regions and the high-contrast regions To address the issue that the match might easily fail over the ice surface with insufficient features, we performed a preliminary operation of distinguishing the lowcontrast regions from the high-contrast regions and set different ranges of template sizes to iteratively determine the optimum template size for each matching position. Experiments were carried out on the single pair of images (WRS2 Path 129 and Row 109) to evaluate the effectiveness of the preliminary operation. Figure 11(a) shows the enhanced imagery and Figure 11(b) shows the result after distinguishing the low-contrast regions from the highcontrast regions, where the low-contrast regions were rendered in red and the high-contrast in blue. As shown in Figure 11(c,d), the match results can be visualized, where the gap areas are the positions without matching. Experimental results suggest that compared with non-image entropy, the proposed method demonstrated a more complete matching results, especially over the low-contrast regions. Some typical areas were selected (as shown in Figure 11(b)) and comparisons of matching completion rate on these typical areas are presented in Figure 11(e). Experimental results indicate that for the high-contrast regions (as shown in Figure 11(b), area 1 and area 2), both the proposed method and the method without using image entropy were able to produce high completion ratio, whereas for the low-contrast regions (as shown in Figure 11(b), area 3 and area 4), the proposed method showed better matching results, with the differences of completion ratio (Jeong, Howat, & Ahn, 2017) over 14%.

Performance comparison with existing data
To evaluate the accuracy of the proposed method, we used the MEaSUREs Annual Antarctic Ice Velocity Maps 2005-2017 (Mouginot & Rignot, 2017) published by the National Snow and Ice Data Center of the United States as reference data. The nominal errors range from 1 m/year along major ice divides with high data stacking to approximately 17 m/year in some areas. We use ArcGIS to register this reference data with the ice shelf velocity field obtained by this paper. The homologous pixels required for registration are mainly located at the stationary areas of the Amery ice shelf, such as, ice cracks, and rock protrusions. In our implementation, 150 checkpoints were selected on the registered images for result comparisons (as shown in Figure 12) and these checkpoints were evenly distributed except the gap areas of the Amery ice shelf. Figure 12(a) presents the ice velocity differences of the selected checkpoints between reference data and that of the proposed method. With the ice velocity values derived from the proposed method as the horizontal axis, and the associated ice velocity values obtained from reference data as the vertical axis, the points were rendered from blue to red according to the value of the velocity difference. Figure 12(b) illustrates ice velocity difference histogram between reference data and that of the proposed method, where the frequency of the ice velocity difference was used as the vertical axis and the associated ice velocity difference as the horizontal axis. We calculated the root mean square error (RMSE) of our method based on reference data and the RMSE of our method is less than 63 m/year. Moreover, the positions with the ice velocity difference less than 5 m/year are presented in Figure 12(c), which suggests that most positions were primarily distributed over the stable areas, such as bare rocks. Consequently, we also concluded that the ice velocity differences between reference data and that of the proposed method appeared probably due to the different time spans and the errors from which reference data might also suffer.

Performance comparison with other existing methods
To further evaluate the robustness and effectiveness of the proposed method, the previous adaptive template size-based method in the spatial domain (Debella-Gilo &  and the globally fixed template size-based method in the frequent domain (Leprince et al., 2007) were compared to the proposed method on image pairs of Landsat 8. For the globally fixed template size-based method in the frequent domain, COSI-Corr based on NCC works within ENVI. Figure 13 illustrates some detailed comparisons between the proposed method and other existing methods, where the gap areas are the positions without matching. Experimental results suggest that the proposed method shows a more complete and more robust matching results than the other methods, especially over the ice surface with insufficient features. Table 1 lists matching results comparisons among different methods, where "total pixels" denotes the number of matching pixels, "matching pixels" denotes the number of right matching pixels and "completion ratio" denotes the matching ratio. Compared with the other methods, the proposed method shows a higher completion ratio with a difference of more than 7.3%. Nevertheless, although the proposed method remarkably improved the matching results particularly over the ice surface with insufficient features, the gap areas still exist (as shown in Figure 13(c,f)) because of the following reasons. First, thin clouds or changing snow conditions might affect the choice of adaptive template size. Second, the ice surface presented significant variants during the time span due to melting snows and other natural environmental changes. Third, the optimum template size might still fail to contain the persistent features between the image pair. As a result, we also confirmed that the performance of the optical image-based matching was directly determined by whether the persistent features on the ice surface existed or not. As mentioned earlier, we performed a multi-scale image matching method based on Gaussian image pyramid for enhancing the matching efficiency. Figure 14 illustrates efficiency comparisons among different methods. For selecting the optimum template size from a set of predefined template sizes, the adaptive template size-based method required a larger amount of computation than the globally fixed template size-based method. Through optimization based on Gaussian image pyramid, the proposed method significantly reduced the computational complexity and improved the matching efficiency.

Discussion
The illustration of the experimental results above is accompanied by brief discussions. Hence, a more general discussion is given here.
In terms of the match able template, we investigate the relationship between NCC and template size. For  the high-contrast regions, the maximum NCC peak can be used for locating the optimum template size by iteratively increasing the template size. However, it is not available for the low-contrast regions since these regions contain less persistent features on the ice surface, which might result in the bad and spurious matching. The low-contrast regions are common on Amery ice shelf, which will lead to many gap areas without the right matching. Thus, the reasonable matching results over the low-contrast regions are of importance for a complete ice velocity field. To distinguish the low-contrast regions and the high-contrast regions, we investigate the relationship between the template size and SNR, and the relationship between the template size and the image entropy. It was observed that the image entropy can be used to distinguish the low-contrast regions from the high-contrast regions when the template size is 23. With regard to the low-contrast regions, the large initial template size that contains more surface features is adopted for finding the adaptive template size, which will increase the probability of a successful match and generate a more robust displacement between the image pair. However, although it can effectively improve the matching results over the low-contrast regions, the proposed method fails to completely eliminate the error matching (as shown in Figure 12(c,f)).
For the image matching method, Debella-Gilo and  pointed out that the computational efficiency needs to be improved due to iteratively finding the locally adaptive template size for each matching position. In this paper, we design a Gaussian pyramid image-based hierarchical data structure to support a coarse-to-fine image matching strategy. Starting from the coarsest image at the top pyramid level, we conduct the improved multi-scale image matching approach and the matching results on the coarseresolution image as the seeds are progressively propagated and expanded to the higher resolution image to obtain a more accurate motion information on the ice surface, which will not only remarkably enhance the matching efficiency, but also effectively reduce the erroneous match positions due to false correlation maxima.
As shown in Figure 15(a), the complete ice velocity field of the Amery ice shelf is derived from the eight pairs of imagery through integration over the overlapping regions and gaps filling. It can be observed that the maximum velocity is mainly concentrated at the inlet of the Amery ice shelf and the intersection of Lambert glacier, Fisher glacier and Mellor glacier. The former is mainly because the ice volume of the three glaciers converges and goes into the sea through the Amery ice shelf, while the latter is due to the topographic relief. In contrast, for the other open and flat regions, the ice velocity is slowed down because of the effect of the wall of the Amery ice shelves. Additionally, it is noted that the ice velocity direction exhibits a great difference between the areas inside Amery ice shelf and those outside the grounding line. More specifically, the direction of ice flow velocity inside the ice shelf exhibits an overall movement trend toward the inlet, reflecting the specific ice flow trajectory, whereas that outside the grounding line does not present such obvious directionality. Some detailed examples of the different positions around the bordering landmass are displayed in Figure 15(b,c), which suggests that the ice flow derived from the proposed method in the vicinity of the bounding landmass also conformed with the overall displacement trend. It further illustrated that the proposed method is robust and effective for retrieving ice velocities from images.
In our future work, the time series analysis based on the ice velocity filed will be used to monitor the change of ice velocity for estimating the mass balance of ice shelves.

Conclusions
In this paper, we proposed an improved multi-scale image matching method for producing a complete and accurate Amery ice shelf velocity field from Landsat 8 images. First, we investigated the relationship between the template size and the image entropy and propose an image entropy-based preliminary operation for distinguishing the high-contrast regions from the low-contrast regions, which could effectively improve the matching results over the regions with the insufficient features. Then, we designed a Gaussian pyramid image-based hierarchical data structure to support a coarse-to-fine image matching strategy and conducted the improved multi-scale image matching approach on the coarsest image at the top pyramid level. During the matching procedure, most computations were performed on the coarsest image at the top pyramid level and the matching results on the coarse-resolution image as the seeds were progressively propagated and expanded to the higher resolution image to obtain a more accurate motion information on the ice surface, which not only remarkably enhanced the matching efficiency, but also effectively reduced the bad and spurious match positions due to false correlation maxima. Finally, we perform a postprocess step for deriving a complete and accurate Amery ice shelf velocity field. Experiments were conducted on the Amery ice shelf using Landsat 8 imagery and the produced velocity field was compared with those acquired from COSI-Corr and MEaSUREs Annual Antarctic Ice Velocity Maps 2016-2017. Experimental results demonstrated that the proposed method can significantly improve the computational efficiency. Moreover, compared with other existing methods, the proposed method can provide more accurate and more robust matching results, particularly over the lowcontrast surfaces. The final results also suggested that the overall trend of the ice flow velocity field generated by the proposed method held good consistency with that acquired from MEaSUREs Annual Antarctic Ice Velocity Maps 2016-2017.

Disclosure statement
No potential conflict of interest was reported by the authors.