A Hierarchical Matching for Amery Ice Shelf Surface Motion Estimation from Repeat Satellite Imagery

In this paper, remote sensing data of Amery ice shelf was used to study Antarctica ice motion and flux problem by a hierarchical image matching method. It combines feature points and grid points to provide a dense, precise and reliable matching result. First, seed points are extracted at the top level of image pyramid using the SIFT algorithm with RANSAC approach to remove mismatches and enhance robustness. These points are used to construct an initial triangulation. Then, feature point and grid point matching are conducted based on the triangle constraint. In the process of hierarchical image matching, the parallaxes from upper levels are transferred to levels beneath with triangle constraint. At last, outliers are detected and removed based on local smooth constraint of parallax. Also, bidirectional image matching method is adopted to verify the matching results and increase the number of matched points. Experiments with Landsat7 images show that the proposed method has the capacity to generate reliable and dense matching results for surface velocity estimation from stereo satellite imagery. Global warming will lead to Amery shelf and glaciers melt and flow rate increase, which can be confirmed by on-site GPS and remote sensing data. Through research the ice shelf flow velocity field, the bottom can calculate the ice flux of this area, and result confirm that the impact of climate for glacier and ice shelf.


INTRODUCTION
The Amery Ice Shelf (AIS) (figure 1-1) is the third largest embayed shelf in Antarctica. It is main thoroughfare of the East Antarctic ice sheet material flows into Atlantic. In the current context of global warming, Amery Ice Shelf's slight changes will exert huge impact on the world environment and global climate. It drains the grounded ice from the interior of the Lambert Glacier drainage basin, which covers 16% of the mass of the East Antarctic ice sheet and is the world's largest glacier by volume [2]. The Lambert Glacier is up to 65 km wide and 400 km long, and drains about 8% of Antarctica's ice sheet to AIS. Today, however, people lack of awareness on the Antarctic ice shelf dynamic deformation characteristics and stability, especially the mechanisms of the surface and inside motion is far from in-depth understanding. Antarctic ice sheet is the past global climate change best record carrier. Antarctic ice sheet with its continuous deposition, the deposition volume per unit time, fewer epigenetic disruption, preservation of the past climate and environmental change resolution, continuous record. This is also the country competed in the Antarctic glacier research carried out mainly. Polar regions to global climate change "Amplifier" role. Historical and modern observations have shown that the magnitude of polar climate change is medium / low latitudes 2 times, indicating that the polar easier to monitor in the middle and low latitudes imperceptible to the subtle changes.
From interferometric analysis of RADARSAT SAR data, Young etc. generated a dense network of surface velocity vectors over the AIS. Ice speeds decrease downstream from about 800 near the confluence of major ice streams at the grounding zone to less than 350 . They then increase to nearly 1400 m a-1 at the centre of the calving front. Young derived the strain rate fields over the AIS from these velocities.
Image matching and feature tracking technique also can be used to measure surface velocity of ice flow. It is to study how to choose some features and similar standards based on the reference image and searching image to search for strategies for correlation computing and to determine the best space responding point for matching. Its main issues focus on the feature space, similarity measurement and searching strategies. The fatal step is to ascertain effective matching methods; a good matching method requires high reliability, small error, fast speed and good real-time.
Image matching is relatively easy when encountered with good image texture conditions. However, on relatively poor textural images, image matching is a difficult and challenging problem. Most of the traditional digital photogrammetry systems require lots of human interactions to remove the errors in the matching results when dealing with poor textural images. SPECIAL FOCUS PAPER A HIERARCHICAL MATCHING FOR AMERY ICE SHELF SURFACE MOTION ESTIMATION FROM REPEAT SATELLITE IMAGERY In this paper, a triangulation-based hierarchical image matching algorithm for stereo satellite imagery is described. It uses a coarse-to-fine hierarchical strategy and combines feature points and grid points to provide a dense, precise and reliable matching result. First, some seed points are extracted at the top level of image pyramid using the SIFT algorithm with RANSAC approach to remove mismatches and enhance robustness. These points are used to construct an initial triangulation. Then, feature point and grid point matching are conducted based on the triangle constraint. In the process of hierarchical image matching, the parallaxes from upper levels are transferred to levels beneath with triangle constraint. At last, outliers are detected and removed based on local smooth constraint of parallax. Also, bidirectional image matching method is adopted to verify the matching results and increase the number of matched points. Experiments with Landsat images show that the proposed method has the capacity to generate reliable and dense matching results from stereo satellite imagery.
In the process of projection from three-dimensional world to two-dimensional images, a lot of information is lost and it is an ill-posed problem. So we must make the best of constraints included in the problem to be solved to limit the size of solution space. After the computation of matching, in reality, there are complex situations such as occlusion, shadows, poor texture, atmospheric dust and steep terrain. So there may be some wrong correspondences and constrains are necessary to eliminate blunders. For now, the common constraints for image matching in digital photogrammetry are mainly as follows: 1) Similarity constraint. Corresponding points are assumed to have similar intensity or colour. So intensity is the main information used in stereo matching.
2) Continuity constraint. Surfaces of objects are assumed to be smooth, which means its parallax varies continuously.
3) Uniqueness constraint. Uniqueness constraint is that searching for matching point on the right image, taking left image as reference, matching point with right image and reference point with left image should be consistent.

II. MATCHING METHODS REVIEW
In the last few decades, a lot of efforts have been devoted in the field of photogrammetry and computer vision to improve the reliability, automation, and efficiency of image matching which can be generally divided into two classes based on the matching primitives. One is area-based matching and the other is feature-based matching.
Area-based Matching: Area-based matching usually works directly on local image windows, and it can acquire dense correspondences. It uses the grey value of the whole image to measure the similarity of two images directly. And a certain method is employed to search the point where the similarity measurement is the biggest. There are many area-based matching methods such as Maximization of Mutual Information, correlation method, conditional entropy method, joint entropy method and so on. Although area-based matching is the most widely used, there exists some shortcomings such as huge computation, long time of matching and sensitivity to rotating, scaling and distort.
Feature-based matching: The common image feature includes point feature, straight line, edge, shape, closed area, statistical moment, etc. By far, feature extraction algorithm can be divided into three main classes: one is point feature extraction operator such as Förstner operator, Harris operator and SUSAN operator, the second is linear feature extraction operator (such as Canny operator, LoG operator) , and the third is surface extraction operator mainly through region segmentation. Generally speaking, feature-based matching has the advantage of being simple to operate, rapid matching speed and high precise matching rate, but it also requires human intervention and the obtaining of feature points is a bit difficult. Besides, it is only suitable for simple images with significant geometric features.
According to the above analyses, feature-based image matching obtains robust but sparse matching results, while area-based matching can obtain dense matching results but the matching reliability may depend on the texture conditions of the images. Therefore, this paper presents a hierarchical image matching method that combines the advantage of both the feature-based matching and the area-based matching methods and produces reliable and dense matching results with high efficiency and automation.
After the matching primitives are selected, the next task is to measure the similarity of the corresponding points in the image pairs, simulating human eyes by means of similarity measurement. The similarity measurement will be the matching score to judge the corresponding points or will be used for a global strategy to judge the corresponding points. There have been many similarity measurement for image matching presented in the past decades, such as the normalized cross-correlation (NCC) (Helava, 1978), sum of squares difference, normalised mutual information (Knops et al., 2006). These are simple computationally, but not distinct enough, and are sensitive to geometrical distortion and discontinuity problems. When using normalized crosscorrelation to measure the similarity between the features on stereo images, an important problem is the selection of an appropriate window size to calculate the correlation values. So in this paper, NCC is selected as the matching score to measure the similarity between corresponding points of a stereo pair and the selection of window size will be discussed in experimental section.

III. DATA AND INSTRUMENTS
A. GPS DATA These data were acquired in the field using a combination of standard surveying techniques (electronic distance measurement and GPS).The velocity vectors are overlaid on a subscene from the USGS 1 km AVHRR mosaic.

B. ICE VELOCITY THROUGH INSAR METHOD
Using RADARSAT SAR imagery obtained during the 2000 Antarctic Mapping Mission, ice velocity vectors were obtained over the Lambert Glacier. The areas of no motion (yellow) are either exposed land or stationary ice. The smaller confluent glaciers have generally low velocities (green, 100-300 meters per year) which gradually increase as they flow down the rapidly changing continental slope into the upper reaches of the faster flowing Lambert Glacier (click on ant-flyover-animation). Most of the Lambert itself has velocities between 400-800 meters per year, with a slight slowing in the middle section. As the glacier extends across the Amery Ice Shelf, velocities increase up to 1000-1200 meters per year as the ice sheet spreads out and thins. Only a handful of in situ velocity measurements have been previously reported of this huge glacier system. While the in situ and radarderived measurements appear to be qualitatively similar, the extent and accuracy of the new measurements are unprecedented and provide quantitative baselines for future comparisons. The ice velocities are obtained from pairs of images obtained 24 days apart, using a technique called radar interferometry. This technique enables a highly precise alignment of image pairs that provides accurate measurements of topography as well as surfaces that have changed or moved over the short time interval, including glaciers.   Table 1: The Landsat7 images used in the analysis GPS data (Table 2) used for validation were acquired in 1991 [8] and 1997-2000 (unpublished) as part of various campaigns conducted over the AIS (see Figure 1-2 for locations).

A. Overview of the approach
The inputs for this approach are the images and rational polynomial coefficient (RPC) parameters. The workflow as shown in Fig. 1   Image pre-processing including image transformation from 16-bit to 8-bit, the Wallis filter and constructing image pyramid.
(2) Feature point extraction using the Förstner. This step is performed to provide the interest points and edges for later image matching.
(3) Coarse image matching using Sift algorithm. This step is employed to generate the initial triangulation for image matching.
(4) Matching propagation in the image pyramid. This step includes feature point and grid point matching.
(5) Blunders elimination in each level including local smooth constraint of parallax and bidirectional image matching.
(6) Least squares image matching in original image level. By means of correcting radiometric and geometric distortion for images, not only is the matching precision improved but also enhance robustness by the rejection of blunder.
This paper will carry on the analysis from the following several important aspects about the approach.

B. Image pre-processing
The quality of matching results relies heavily on the quality of the images. Poor contrast would reduce the number of interest points extracted. The Wallis filter is employed to strongly enhance and sharpen the already existing texture patterns and increase the signal-to-noise ratio, which is beneficial for feature extraction.
At present, the grey values of panchromatic remote sensing images are stored by 16-bit. However, in our research, 8-bit stored images are necessary for the computation of image matching and the display of images in order to decrease calculation amount. So the 16-bit images are mapped to 8-bit images using a principle of auto levels.
After that, we must construct the image pyramid before matching. There are two kinds of methods for generating image pyramid: the direct method and the filtering method. The direct method is to transform 2 pixels by 2 pixels to 1 pixel with an average value. Filtering method, which is adopted in this paper, varies from the first one in that low-pass filtering is used to replace the average. Specifically, starting with the original image, each subsequent level of the image pyramid is created by sub-sampling the previous level image and smoothed by a Gaussian filter. Lower level images have coarser resolution and details are lost due to smoothing. The pyramid level number is a pre-defined value which could be either a user-input or can be determined according to the height range of the imaging area.

C. hierarchical image matching
The hierarchical searching algorithm is by narrowing the search scope to achieve the purpose of reducing the computational complexity to improve matching speed. It is proposed just as people find things with the coarse-tofine strategy.
The hierarchical image matching method firstly uses a SIFT algorithm and RANSAC approach to obtain a few reliable correspondences, and then construct an initial triangulation. The SIFT algorithm (Lowe, 1999) is proved to be able to produce robust but relative sparse corresponding points invariant to moderate scale changes or distortions, which is ideal for the purpose of generating a certain number of well distributed matching points for the initial triangulation. In the SIFT descriptor, each interest point is characterized by a vector with 128 unsigned eight-bit numbers generated from a local region, which defines the multi-scale gradient orientation histogram. The matching is performed by measuring the similarity between the two vectors associated with the two matching points.
The RANSAC approach is used to detect and eliminate possible blunders from the previous SIFT matching results. It starts by randomly selecting a subset of the matched corresponding points. From the chosen matched points, a fundamental matrix can be calculated based on which a model is then built. This model is evaluated by determining whether each pair of corresponding points fit reasonably well to it. This is used as a criterion to determine the best model which has the largest number of correct corresponding points. This process is repeated to find the best model. Those matched points which do not fit for the final best model are considered as blunders and eliminated from the initial matching point set (Wu, et al., 2012).
After the seed points are extracted at the top level, an initial triangulation can be constructed and an area-based image matching with feature points and grid points is conducted at the top level again.
In the process of hierarchical strategy, image matching is first conducted on the lowest resolution. The matched points are then transferred to the next level (of higher resolution) where additional feature points could be matched. This process repeats until it reaches up to the original image level. At a subsequent level, points from upper level are matched again to achieve higher precision. A TIN (Triangulated Irregular Network) surface of parallaxes is generated from these matched A HIERARCHICAL MATCHING FOR AMERY ICE SHELF SURFACE MOTION ESTIMATION FROM REPEAT SATELLITE IMAGERY points using the Delaunay triangulation. This TIN is used to estimate the correspondence of additional feature points (Hwangbo J. W., 2010). Interpolation of x parallax using TIN surface As shown in Fig. 4

D. Feature point and grid point matching
Förstner operator is one of the most famous operators to extract feature point and it has high speed and positioning accuracy. Förstner operator generates a series of feature points by identifying distinctive variations of pixel brightness values. The detected feature points including corners and round dots can reflect the shape of the terrain in the images. Feature points are typically found around shadow, where significant change of brightness values occurs. These points are related to the shape of terrain or objects such as ridges and rocks.
Feature point matching is very efficient and suitable in texture-rich regions with grey value variation. On the other hand, in image regions with poor texture or no texture information, few or even no feature points can be extracted. Thus, feature point matching will lead to holes in these areas. To solve this problem, grid points can be used and grid point matching has been introduced.
Grid points are points determined at given positions that are generally uniformly distributed over the whole image. Compared to the feature points, the choice of grid points is blind and thus many grid points may lie in poorly textured regions or occluded areas. The search for the match of a given grid point has a higher possibility of yielding an ambiguous match and even no matching candidates. The grid points are also matched using crosscorrelation following the coarse-to-fine strategy.

E. Blunder elimination
Even though the constraint of epipolar line and triangle is adopted in the process of image matching, still it cannot guarantee all the matching results are right. So it is important to employ an algorithm for detecting and removing outliers.
At first, bidirectional image matching algorithm is imbedded in our program. The basic idea of this strategy is simple as follows: First, a point in left image is taken as target point and find out the point of maximum similarity in right image. Then take the matched point (the point in right image) as target point and search for the best point in left image in the same way. Finally, we accept it if the two results overlap, or give it up as mismatch.
Although bidirectional image matching method eliminates some wrong matches, there still exist some mismatched points. Therefore, local smooth constraint of parallax is adopted in this research. The logic behind the strategy is that the terrain is continuous and spatially correlated in a small area. Speaking specifically, for every pair of corresponding image point, the adjacent corresponding points within a certain distance are used to fit an optimal plane of parallax. Outliers are the values that are statistically far from most others in a set of data. Since the spatial distribution pattern of the parallax of corresponding points could vary depending on the type of local terrain, the parallax offset from the fitted plane is compared with the standard deviation of the neighboring points. If the offset is beyond 3 , the matched point is considered as a blunder and discarded from the result.
V. EXPERIMENT RESULT We have applied the proposed matching approach to Landsat stereo images for ice surface motion measurement in AIS, , East Antarctica. The resolution of the images is 15 m. The scene covers a total area of 175km by 350 km and consists of a variety of land cover types, including blue ice, ice-free rock, crevasse, snow and firns, water body and moraine. Controlling the matching parameter setting is important for high quality matching generation. Since the stereo matching result is highly dependent on the properties of input data, there is no single set of parameters that is perfect for every image (Hwangbo J. W., 2010). So Table 1 shows some comparatively optimal parameters we selected through many experiments.  They provide robust structure of terrain from coarse level with a few points to finer level with an increased amount of detail. Fig. 5-3 illustrates the generated HMV(Hierarchical Matching Velocity) with Rignot published velocity field. Because many areas have cloud cover or significant changes, there are many blunders in these areas. Figure 5-4 illustrated the velocity field. This velocity field is generated from the matched points using triangular linear interpolation approach. Visual inspection reveals that good results have been achieved. selected and checked to ensure that these automatically matched velocity. In this area, points were manually matched in the stereo images for generating a manual matching velocity. Then GPS, MEaSURE and manual matching velocity was compared with the HMV generated from automatically matched points in the of profile.. Overall, the velocity result from 4 different methods showed high consistence. The differences between manually and automatically generated velocity field were in the range of 25 maximum with standard deviation less than 6m. The random elevation errors of such magnitude are reasonable, considering that the pixel resolution of the raw images is about 15m. vertical direction of the velocity and flow of various points can be approximated as identical, and secondly, the flow rate of each point are essentially the same, are similar to the cross-sectional direction; once again, at the bottom did not like the undulating ground ice, which can be considered as approximate plane. Based on the above conditions, we can establish the model shown in Figure  6-1. Therefore, the calculation section ice flux is actually calculated by the cylinder as shown in irregular quality of the ice contained, if the ice flow and cross flow direction is not perpendicular to the section, the flux result need to multiply sin!. Section i of the ice flux is calculated as " " In the formula, ! is the angle of the ice flow direction between the section.. Ice at the bottom of the ice surface velocity of flow is 0.8 times, and the ice surface to the bed of ice flow on the vertical variation of uniform size, equivalent to the entire surface of the cylinder, the average velocity of the flow rate of 0.9 times. Figure 6-1 Therefore, the calculation of the Lambert Glacier Basin, ground line output flux material requirements for the major InSAR ice velocity data, the ice thickness data, ice density data. Ice thickness data is the use of hydrostatic fluid balance formula, derived from the Amery Ice Shelf DEM obtained. Because the data distribution and spatial interpolation error is generated when the DEM, geoid model OSU91a (error of about 3m) and ice density data uncertainties (error of about 5kgm-3), calculated by the DEM in some parts of the ice thickness 100 ~ 200 m possible errors. The flow rates of data error of about 5-10 m a-1. Taking all the above errors ice flux across the grounding line total error is set to 10% (Wen et al., 2007). It can be concluded that the entire Lambert Glacier Basin, ground line ice flux of about 88.9 ± 8.9 Gt/a.

VII. DISCUSSION AND CONCLUSIONS
This paper presented a triangulation based hierarchical image matching method for stereo satellite imagery. Table 2 gives the accuracy evaluation results in the whole rectangular region mentioned above. About 94% of the region is within three standard deviations though the maximum absolute difference is 4.841m.
Though the experiment using Landsat7 ETM images, we draw the following conclusions: (1) The SIFT algorithm with RANSAC approach can produce reliable but sparse corresponding points which is ideal for obtaining the initial values for hierarchical image matching.
(2) The hierarchical matching strategy based on triangulation constraint proved to be effective and reliable especially in some area with poor and repetitive texture. More importantly, for stereo satellite imagery, the hierarchical search strategy can improve the speed of image matching; otherwise it will be very time consuming to search corresponding points in the original image.
(3) The experiment results showed that local smooth constraint of parallax can eliminate most of outliers effectively, although there still exist some mismatches. Future work will be focused on this problem to better refine the matching results.
The widespread retreat of glaciers can be considered as a response to the climate change. Being the largest retreating glacier-ice shelf system in East Antarctica, the Amery Ice Shelf-Lambert Glacier system plays an important role in contributing to sea level rise as well as the surrounding environment and climate. The present study is focused on the investigation of the stability of Amery ice shelf and Lambert glacier system, mainly surface ice flow field flux calculated. Global warming leading to increased Antarctic glacier change, flow speed up, greater flux. Article does not consider the surface melting and freezing, because studies show that the main mass lose through under ice melt and glaciers flow. Future, if considering the cumulative rates of sedimentation can get estimates of the mass balance .
ACKNOWLEDGMENT I would like to gratefully acknowledge the enthusiastic supervision of Dr Ron Li during this work. I thank Prof Weian Wang for the technical discussions. This document's development has been supported by the polar 973 project, and peer reviewed by Associate Professor Qiao.