A NOVEL METHOD FOR ESTIMATION OF GLACIER SURFACE MOTION IN 1960 s FROM ARGON KH-5 OPTICAL IMAGERY

The Antarctic ice sheet response to the global climate change, specifically the ice flow speed change of the glaciers, has been investigated by many researchers. However, most research results cover the period since 1970s or after the operation of the LANDSAT series. The availability of the film-based ARGON KH-5 data makes it possible to quantify the changes of the Antarctic ice sheet in 1960s. To meet the challenges of processing the low quality film-based ARGON images, a novel method was developed to allow estimating the ice sheet surface motion and reconstructing the surface model simultaneously from ARGON stereo images by decomposing the total parallaxes to terrain and motion based components. A photogrammetric approach was developed to distinguish stable ice surface features from those on motion and use them for recovering the camera orientation information. Several existing Antarctic mapping products were used to establish the ground control. The ice flow speed field is reconstructed using a hierarchical image matching strategy. Firstly, epipolar images are generated via a fundamental matrix derived from correspondences used in the geometric modelling process, and then an image pyramid is built. Second, the normalized cross-correlation (NCC) technique is conducted on each layer of the pyramid to match the extracted features. Since the images were taken at different times, during which the glacier motion occurred, the measured total parallaxes are decomposed to terrain and motion parallaxes according to given ice flow directions which are derived from the iteratively produced DTM or images. Finally, a speed map and a DTM can be generated at each level of the image pyramid. This process repeats itself. At the bottom of the pyramid the final speed map and DTM are produced at a resolution of about 60m and represent the ice flow field of 1963. This approach was tested using two ARGON stereo-pairs in Rayner glacier in East Antarctica. Both the ice flow speed map and DTM were generated, and their difference with recent products is briefly discussed. * Corresponding author


INTRODUCTION
Effects of global climate change have been repeatedly reported for the past decades.As one of the significant responses to the global change, the ice velocity has been investigated to estimate the Antarctic ice mass as well as its contribution to the global sea level change.In-situ measurement is a direct way to obtain the ice flow speed information.Various approaches, such as stake measurements (Budd et al., 1982) and GNSS observations (King et al., 2007;Urbini et al., 2008), were applied to derive the glacier surface motion.Yet it is not easy to accomplish the ice velocity estimation for a large area given the time and effort it would take to perform the observation.Alternatively, the remote sensing-based estimation using radar interferometry, optical images and laser altimetry (Lee et al., 2012;Liu et al., 2012;Hulbe et al., 2013) serves as an excellent complement.Compared with the ground-based method, the remote sensing approach can produce timely observing results over regional scale that is not reachable by human beings due to the harsh environment, making it a major means in research field to obtain the Antarctica ice flow velocity during the recent decades with modern satellite observations.Early velocity measurements of Antarctica ice sheet, which is derived from historical remote sensing records, play a vital role for estimating ice flow velocity variations and reducing the possible fluctuations over long time scales.However, most research results only cover the period since 1970s or after the operation of the LANDSAT series.The availability of the film-based aerial photographs and declassified intelligence satellite imagery (DISP) including CORONA, ARGON, and LANYARD missions (Ruffner, 1995) makes it possible to quantify the changes of the Antarctic ice sheet in 1960s.All the films were retrieved via re-entry capsules, which were separated from the satellite and fell to earth.Figure 1 illustrates the recovery maneuver of films.The previous work related to the early film-based photographs mainly focuses on mathematical modelling using DISP data (Zhou et al., 2002;Kim, 2004;Sohn et al., 2004).Related issues yielded during the geometry reconstruction, such as film image scanning (Gheyle et al., 2011), image pre-processing (Galiatsatos et al., 2005) and DEM generation (Surazakov and Aizen, 2010), were investigated as well.The research for employing historical aerial and satellite imagery in ice sheet studies emphasizes the surface change analysis with the assist of modern images.For example, historical aerial and satellite photographs were processed to depict the retreating trends in 224 margin glacier fronts in the Antarctica Peninsula and associated islands ever since 1940s, revealing patterns broadly compatible with retreat driven by atmospheric warming and other drivers, and also an overall reduction in total ice-shelf area of 28,000 km2 since the beginning of the period was found (Cook et al., 2005;Cook and Vaughan, 2010).Furthermore, aerial stereo photographs acquired in 1960s were processed to extract DEMs and then subsequently compared to DEMs derived from modern ASTER satellite data and aerial photography.Volume changes for two glaciers in Antarctica Peninsula were measured and analysed (Kunz et al., 2012).
Currently, by reviewing the existing literature, there is no Antarctica ice velocity product available that derived from the film-based ARGON DISP images in 1960s.This attributes mainly to the challenges for processing the declassified satellite photographs.This paper is to propose a new analytical method for estimation of Antarctic ice sheet surface motion in 1960s using historical DISP photographs.This method is based on the parallax decomposition of two adjacent matched ice-flow features from two DISP images with different acquisition time.During this process, parallax constraints were added as auxiliary information to help matching.Finally, besides the ice flow velocity, the terrain model on the ice flow region could also be derived from the parallax analysis as a by-product.

ICE FLOW VELOCITY ESTIMATION FROM SEQUENTIAL PERSPECTIVE IMAGES
A novel method is developed to allow simultaneous estimation of the ice sheet surface motion and surface terrain model from merely one ARGON stereo images acquired at different times.
The idea of the proposed ice flow estimation method is based on parallax decomposition of two adjacent matched ice-flow features from two images acquired at different times.The matched features are derived from the image matching process in which a hierarchical matching strategy is involved.The details of the method are explained below.

Ice Velocity Estimation via Parallax Decomposition:
This approach is feasible if and only if two essential conditions are met, that is, the stereo image planes should be almost parallel to the ground plane, and the stereo pair from which the correspondences are identified should be transformed to epipolar images.Both of them significantly simplify the geometry relationship between the ground space and the image space.Once they were satisfied, the parallax of matched iceflow features would only be affected by the terrain and the surface motion.Specifically, two components are contained in the horizontal parallax of ice-flow matched features, i.e., the one resulted from the terrain and the other one resulted from the surface motion.As for the vertical parallax of matched features, it should be close to zero theoretically if the features are static, yet a significant shift in vertical parallax can be observed if they moved.In other words, the vertical parallax of matched ice-flow features can be considered as the vertical component of the motion in image space.Thus, if the motion direction can be estimated, the horizontal component of the motion in the horizontal parallax can be separated.Actually the motion direction can be accessed through the aspect of the terrain.As a result, the total ice motion can be recovered from the decomposed parallax as well as the terrain information.

Image Matching:
Matching is the key step to the ice velocity estimation, and the accuracy of the ice velocity is up to the accuracy of matched correspondences to a great extent.The overall matching steps goes as follows: Firstly, some accurate, well-distributed correspondences, which were selected as tie points in bundle adjustment process, are used to get the fundamental matrix through random-sample consensus algorithm (more than 8 corresponding points are needed) to generate epipolar images (Richard, et al, 2003).Moreover, these points are taken as check points to verify the accuracy of epipolar images.Subsequently, Gaussian image pyramids are established for both stereo epipolar images and image matching is carried out through the top to the bottom level of the image pyramid based on the normalized cross correlation (NCC) to obtain correspondences.In the process, triangulated irregular networks (TINs) generated from matching control points are employed to provide constraints for identifying both the horizontal and vertical parallax of potential correspondences.Newly matched features then are added into the control point set to update the constraints.Finally, this procedure repeats itself till the last level of the image pyramid and a dense dataset of correspondences could be expected.It is note worthy that, to avoid introducing extra errors caused by the correction for distortions, e.g.image warping, the matching process is carried out first using images without corrections.Once the correspondences were identified, then triangulation is performed with corrections to obtain the precise positions of these correspondences in the ground space.

Area and Data
The study area is Rayner Glacier, a prominent glacier about 16 km wide, drawing part of Enderby Land in East Antarctica, flowing north to Casey Bay, with ice flow speed above 100 m/year for more than 200 km inland, and the velocity can even speed up to 1,000 m/year on the Rayner/Thyer Ice Shelf (Rignot, 2011).The recent overall mass change of the ice shelf was almost balanced, but a relatively significant basal melting was reported (Liu et al., 2015).
Two stereo DISP images covering the Rayner Glacier and Rayner/Thyer Ice Shelf are used in our study for ice flow velocity estimation.They were obtained from two missions, i.e., 9058A and 9059A, which carried the similar camera with a focal length of 3 inches (approx.0.0762m) and were operated in August and October 1963, respectively.Each photograph had the similar format size of 4.5 inches×4.5inches (approx.11.43cm×11.43cm),with a film resolution of 30 line-pair/mm and a ground resolution of approximate 140 m (Light, 1993).
The ground control points (GCPs) were selected from existing Antarctic digital products, namely, the 30m resolution Landsat Image Mosaic of Antarctica (LIMA) and the 30m resolution ASTER Global Digital Elevation Model (ASTER GDEM).The acquired GCPs were further confirmed by visual check of the red-cyan anaglyph 3D image and a within-grid selection strategy was applied to ensure their reasonable distribution among the images (Ye et al., 2016).

Pre-processing
In view of the special characteristics of the ARGON DISP images and the overall process of the proposed method, the preprocessing section is composed of two parts, one is image denoising and enhancement for better feature extraction and matching, and the other is orientation parameter determination through bundle adjustment to recover the rigid geometric relationship.
The noises associated with the ARGON DISP images have different types, including the random speckle noises caused by low-quality exposure when imaging, aperiodic stripes and scratches introduced by long-term reservation, and Newton's rings generated during the scanning process.Several de-nosing methods were tested to eliminate these kinds of specific noises and some experimental results were obtained (Ye et al., 2016).
The second part focuses on the camera pose determination.An essential step is to make approximations for parameters of both the interior orientation (IO) and exterior orientation (EO).The IO elements can be derived from measured positions of fiducial marks in the image and the corresponding calibrated coordinates derived from the camera calibration report.The 2nd order polynomial model was adopted in the process to depict the transformation relationship.for the EO elements, they can be estimated via the Lagrange polynomial using the satellite ephemeris and the ground control.

Results and Discussions
After estimating the camera poses, bundle adjustment provided by the Leica Photogrammetry Suite software was performed with 23 tie points and 5 ground control points, where all the points involved were located either on distinctive bare rocks with no motion or slow moving features with a speed lower than a 100 meter per year.The total image unit-weight root mean square error (RMSE) of bundle adjustment turns out an accuracy of nearly half a pixel, indicating a good internal accuracy of the adjusted model.Besides, 5 check points were selected in the same criteria as GCP to examine the external accuracy of the bundle adjustment, and it showed an average ground residual under half a nominal pixel (140 m).
The tie points and GCPs used in the geometric processing were also used to generate the epipolar images.Another thirty distinctive stable features were selected to verify the quality of the epipolar images and the mean vertical error is under one pixel.The constrained NCC matching was performed through the image pyramid.The terrain model of the ice sheet was produced on each layer.However, the parallax decomposition operation was conducted from the third layer given the significance of the ice motion influence on the vertical parallax.Thus, only four ice flow speed field maps were produced.By comparison of either the terrain model or the ice velocity map from different layers, it can be found that as the increase of the correspondences, the pattern of both products turned smooth and obvious.
By comparing this speed map with a recent one by Rignot (2011), it shows that the overall speed in 1963 was higher, though the spatial speed pattern appears to be similar, including both the main glacier and tributary ice flows.A significant advance of the ice shelf front can be observed after georeferencing.Besides, the obtained surface model was compared with ASTER GDEM after motion correction, and it shows a mean elevation difference of only a few meters.Further detailed analysis of the speed changes vs. the Rayner environment is being conducted.