DEM RECONSTRUCTION OF LARSEN B REGION BASED ON 1960 S OPTICAL SATELLITE IMAGERY

Antarctica is an important part of the earth system and crucial to global sea-level and climate change. In 1995, a number of US intelligence satellite photographs from the1960s to 1970s have been declassified, only a portion of which (ARGON) covers the Antarctica. It provides a broader perspective to study the early Antarctica and it’s very important for us to study the changes in the early stage of Antarctica. In this paper, a hierarchical stereo image matching strategy was used to reconstruct a digital elevation model of the Larsen B region of Antarctica Peninsula in 1960s using ARGON images. The accuracy of the matching result estimated in all layers is within one scanned-pixel of 33m and the accuracy of geometric modelling after bundle adjustment estimated by using check points is within one nominal pixel of 140m. In the future, the elevation changes from 1960s of Larsen B tributary glaciers will be analysed.


INTRODUCTION
Antarctic Peninsula (AP) has undergone significant changes in the area of ice shelves, ice thickness, air temperature and et al. over the past few decades.Studies on AP shelf changes have shown that the regional environment changes can cause substantial losses in ice shelf (Cook and Vaughan, 2010).A lot of researches have been published about Antarctic Peninsula, especially Larsen B ice shelf.With the rise of regional climate in the past decades, accelerated retreat and break-up have been observed on Antarctic Peninsula ice shelf (Vaughan and Doake 1996, Skvarca et al., 1999, Rack and Rott, 2004).Skvarca et al. (1999) analysed 34 year satellite time series data from 1963 to 1995 suggested Larsen B lost 2320 km 2 in area at the end of January 1995, Rack and Rott (2004) summerized the area of Larsen B ice shelf changes between 1995 and 2003, during that period the ice shelf area decreased from 11512 km 2 to 2667 km 2 , only 2400 km 2 remaining in 2009 (Cook and Vaughan, 2010).On eastern Antarctic Peninsula occurred two big collapses, Larsen A ice shelf in 1995 and Larsen B ice shelf in 2002.A new NASA study reveals that remnant Larsen B ice shelf likely to disintegrate completely by the end of the decade (Steve Cole and Alan Buis, 2015).Besides the area change, numerous studies focused on ice velocity and ice thickness changes are published.The largest tributaries of Scar Inlet ice shelf, Flask and Leppard glaciers, decrease rate 38% and 46% higher in 2013 than 1995 (Wuite et al., 2015).From 2000 to 2003, Hektoria, Green and Evans glaciers increased eightfold and declined since 2003; Jorum and Crane glaciers accelerated continuously from twofold in 2002 to threefold in 2003; on the contrast, farther south Flask and Leppard glaciers had no significant changes (Rignot et al., 2004).Khazendar et al. (2014) used RANDSAT-1 (1997RANDSAT-1 ( , 2000RANDSAT-1 ( , 2004)), ALOS PALSAR (2006PALSAR ( -2010) ) and TanDEM-X (2012) to obtain the velocity of Scar Inlet ice shelf.The result showed the ice flow speed accelerated from 475 m a -1 to 750 m a -1 during 2003and 2013. From 1997to 2012, Flask and Leppard glaciers increased by 55% and increased from 325 m a -1 in 1997 to 400~500 m a -1 in 2004.From 2006 to 2012, the acceleration rate reduced and speed was about 450 m a -1 to 520 m a -1 .Shepherd et al. (2003) revealed that averaged surface decline trends of Larsen B ice shelf is -0.17±0.11m a -1 during 1992 and 2001.After calculating the elevation change rates of six major glaciers of Antarctic Peninsula, Fricker and Laurie (2012) concluded that this area was lowering in average rate of -0.03~-0.16m a -1 thirty years after 1978 and average dh/dt was about -0.2 m a -1 in six inter-mission crossovers from 1992 to 1997, two of which keep the rate of ~-0.19 m a -1 during 1998 to 2008.By employing the fluxes difference between the pre-and post-collapse, Rott et al. (2011) estimated the mass loss of Larsen B embayment, for example 4.34±1.64Gt a -1 in 2008.The combined mean loss rate of Larsen B tributary glaciers during 2001 to 2006 was estimated at least 11.2 Gt a -1 by tracking ice extent and elevation changes (Shuman et al, 2011).Berthier et al. (2012) evaluated the mass loss rate of Larsen B tributary glaciers was (9.0 ±2.1 Gt a -1 ) during 2006-2010, closing to the rate (8.8 ±1.6 Gt a -1 ) from March 2002 to 2006, derived using DEM differencing and laser altimetry.
There is very little data to estimate ice thickness and surface elevation changes before 1970s.The 1960s ARGON optical satellite images declassified in 1995 provide valuable historical data for studying the early changes in the Antarctica.In this paper, one pair of ARGON stereo images in the same track were used to reconstruct three-dimensional terrain of Larsen B region by hierarchical matching strategy.The study is important for revealing the elevation changes and mass balance of the area over half a century.

METHODS
This section focuses on the process of DEM establishing, including the epipolar images generation, the matching process and the bundle adjustment.The ARGON images used to establish the Larsen B DEM are DS09058A014MC115 and DS09058A014MC115.

Epipolar Images Generation
Matching based on epipolar images could reduce the search ranges in the following stereo image matching process.Since Helava (1972) proposed the concept of epipolar in 1970s, a lot of ways to establish the epipolar image were emerged.Here we used fundamental matrix calculated from a series of corresponding image points on two images to perform epipolar images.The steps of the epipolar images generation for the stereo image were discussed in details by Hartley and Zisserman (2003).

The Matching Process
Stereo image matching is one of the key techniques to establish the DEM.A novel method performed to estimation of terrain information and motion speeds from a stereo image pair simultaneous was proposed (Li et al. 2017), the matching strategy in this paper referred to this top-down hierarchical matching method.The study area is 6046×5458 pixel size at original image layer.In order to avoid the effects of tiny features and image noise, a five-layer Gaussian image pyramid was applied to epipolar images.The details of matching process were as follows.

1)
In the 1st layer, SIFT (Lowe 1999) operator was used for feature detection and matching.In generally, the distribution of matching points is critical to sketching the outline of topography in the entire study area, so some points on the top of the mountains and the foot of the mountains was added manually (Figure 1).
2) The SIFT matching points and manual points were derived from step one as seed points (Row 1 and first layer in Table 1) to generate a triangulated irregular network in both images respectively, and then using Shi-Tomasi (Shi 1994) corner detector to detect the feature points (Row 2 in Table 1) in both images.According to the position of the feature point in the left triangulation, the position in the right triangulation is predicted according to the parallax of the three vertices of the left triangulation.Feature points within a certain window around predicted position on right image were regarded as potential matched points, whereas the same.After feature-to-feature normalized correlation coefficient (NCC) matching under the constraints of the triangulated irregular network was performed, a reasonable NCC threshold was used to remove the mismatched points (Row 3 in Table 1).As a result, the seed points and the matching feature points (Row 6 in Table 1) were the totally corresponding image points in 1st layer (Row 7 in Table 1).3) The matching points were passed from upper layer to next Layer (Row 1 and second to the fifth layer in Table 1).According to the down sampling interval when the pyramid image was created, that one pixel on the higher layer was divided into four pixels on the next layer.So the coordinates of the matching point should multiply two when passed the matching points from upper layer to next layer.When one match point transmitted to the next level, four points that divided from passed right matching point were potential corresponding points of passed left matching points, so it is necessary to find the unique corresponding points of the left match point by re-matching.
4) The refined matching points can be as the seed points in this layer to generate a triangulated irregular network in both images respectively.Repeat step 2 and 3 until the match points were passed to the original resolution image.The matching points in each layer and eliminated points by NCC threshold, parallax constraint (Row 4 in Table 1) and surface fitting constraint (Row 5 in Table 1) are shown in Table 1.

Bundle Adjustment
Figure 2. The location of Larsen B with an ARGON mosaic as background (Kim et al. 2007).And the distribution of tie points, ground control points and check points In order to project the matching point to the object coordinate system, it is necessary to calculate the exterior orientation element of images.In this paper, Bundle adjustment (BA) was performed to calculate the exterior orientation element of ARGON images using the commercial software (ERDAS, 2013) LPS of Leica Systems.In this paper we use two ARGON stereo images covering the Larsen B ice shelf.Figure 2 shows the coverage of two ARGON images.22 tie points, 6 ground control points (GCPs) were involved in the bundle adjustment process (Figure 2).The GCPs obtained from LIMA and ASTER DEM, in order to get a good geometry to avoid unpredictable distortion, it is necessary to make sure the distribution of the tie points is balance in the overlapping area.The internal and external accuracy of BA results was evaluated using the tie points and check points, respectively.The result show that the bundle adjustment of the Larsen B yielded a unit-weight RMSE of 0.52 pixels, this indicated a high level of internal accuracy of the BA.
With the geometric positioning parameters obtained from the BA, 5 check points were calculated to assess the external accuracy of the BA.The check points were calculated from BA is compare to the coordinate measured in the LIMA mosaic (Bindschadler, R. et al., 2008) and the ASTER GDEM (Tachikawa, T. et al., 2011).The differences between BA and known coordinates were used to evaluate external accuracy shown in

Larsen B DEM Reconstruction
The matching points were projected onto the world coordinate system using the parameters obtained from the BA.According to the density of the point a 300 m resolution DEM was produced.A 3D view of the ARGON image draped on the DEM shows in Figure 3.

Figure 1 .
Figure 1.Distribution of SIFT feature points and manual points; red dot indicate SIFT feature points, white dot indicate manual points

Figure 3 .
Figure 3. ARGON image draped on the DEM3.EVALUATION OF MATCHED POINTS

Table 1 .
Matching and eliminated points for each layers in details Finally, dense matching in original resolution image was applied to very five-pixel grid.

Table 2 .
We can see both horizontal accuracy and vertical accuracy are within one nominal pixel of 140 m.

Table 2 .
The external accuracy of Larsen B BA estimated by check points