Innovative Method of Combing Multidecade Remote Sensing Data for Detecting Precollapse Elevation Changes of Glaciers in the Larsen B Region, Antarctica

The Antarctic Peninsula has undergone dramatic changes in recent decades, including ice-shelf melting, disintegration, and retreat of the grounding line. The Larsen B ice shelf is of particular concern due to the unprecedented ice-shelf collapse in 2002. Since few observations on the Antarctic Peninsula were available before the 1970s, long-term investigation of the surface elevation change in the Larsen B region could not be pursued. In 1995, the United States administration declassified a collection of archived intelligence satellite photographs from the 1960s to the 1970s, including analogue satellite images from the ARGON program covering parts of the Larsen B region. We chose overlapping ARGON photos captured in the Larsen B region in 1963. These photos were all subjected to a tailored photogrammetric stereo-matching process, which overcomes those specific challenges related to the use of historical satellite images, such as poor image quality, low resolution, and a lack of high-precision validation data. We discovered that between 1963 and 2001, the surface elevations of the main tributary glaciers in the Larsen B embayment have undergone little change before the ice shelf collapse from 1963 to 2001 by comparing the reconstructed ARGON-derived digital elevation model (DEM) (1963) and ASTER-derived DEM (2001). In addition, the results demonstrated that the hierarchical image matching method can be modified and applied to reconstruct a historical Antarctic DEM using satellite images acquired ∼60 years ago through an innovative and rigorous ground control point selection procedure that guarantees no changes occurred at these points over the period. The new ARGON-derived DEM derived from ARGON (1963) can be used to build a long-term spatiotemporal record of observations for extended analyses of ice-surface dynamics and mass balance in the Larsen B region.


I. INTRODUCTION
A CCORDING to the current research on the Antarctic Peninsula, changes in climate and ocean conditions may have resulted in significant mass losses in this region [1], [2]. Because of increased atmospheric warming and ice-shelf basal melting in recent decades, several ice shelves surrounding the Antarctic Peninsula have shrunk and broken apart [3], [4], [5], [6]. Changes to the Larsen B ice shelf (LBIS) and tributary glaciers (see Fig. 1), such as glacier thinning, ice-flow acceleration, and ice-shelf collapse, have resulted in increased mass losses and contributed to accelerated sea-level rise. These changes indicate regional instability, which also has an impact on the stability of the cryosphere and the global climate system.
The LBIS experienced two catastrophic disintegration events in 1995 and 2002. It lost 2320 km 2 of ice at the end of January 1995, after a relatively stable period between 1963 and the warm summer of 1995 [4], [6]. From January 1995 to early February 2002, the LBIS experienced a series of sustained and stable retreats, followed by a rapid and dramatic collapse in the period from February to March 2002, with a total areal loss of 5729 km 2 [6]. These phenomena were caused by an increase in atmospheric temperatures during the summer of 2001 and 2002 [7]. Cook and Vaughan [1] estimated that only 2400 km 2 of the LBIS remained in 2009, compared to approximately 12000 km 2 in 1963.
The disintegration of ice shelves has the potential to cause immediate accelerations in upstream glaciers' movement and thinning [8], [9], [10], [11]. Using ERS radar altimeter measurements, Shepherd et al. [5] calculated an average decreasing surface elevation trend of −0.17±0.11 m/yr in the LBIS between 1992-2001. There were no systematic elevation observations on a decadal scale made over the LBIS tributary glaciers before the 2002 collapse from 1963 to 2001. The ice-flow velocity of the Hektoria, Green, and Crane glaciers increased between 2000 and 2003, with the Hektoria glacier experiencing a six-fold increase in speed and a 38-m drop in surface elevation in six months after collapse [10]. Following that, Pritchard et al. [12]  The southern glaciers, on the other hand, did not change significantly in velocity and elevation. As a result, it was assumed that the Scar Inlet Ice Shelf, a remnant of the LBIS, supported these glaciers [10]. Previous study estimated the elevation changes of 12 glaciers in another region (northwest part) of the Antarctic Peninsula [14]. From the mid-1960s to the 2010s, an annual mean lowering rate of 0.28±0.03 m/yr was discovered using aerial stereo photogrammetry from airborne photographs [14].
Following the 2002 collapse event, due to ice shelf front calving, mass loss increased significantly, resulting in a negative mass balance state in the Larsen B region. From 2002 to 2006, the mass loss rate for the Larsen B tributary glaciers was 8.8±1.6 Gt/yr on a near-decadal scale and 9.0±2.1 Gt/yr from 2006 to 2010 [15]. The total mass deficit of the Hektoria, Green, Evans, Jorum, and Crane glaciers was 21.4 Gt/yr in 2003 [9], which then decreased to 4.3 Gt/yr from 2008 to 2009 [16].
Surface changes will reflect sustained mass loss over a relatively long-time period. The observed surface lowering rates on a decadal scale in the northern Antarctic Peninsula were faster than the ones in the southern Antarctic Peninsula [14]. The processes driving these changes did not feature a simple genesis due to the high local variability [17]. Since the 1990s, within the Larsen B drainage basins, extensive break-up and subsequent changes were observed [4], [6]. The comparison of multitemporal digital elevation models (DEMs) allows for quantifying surface elevation changes. DEMs in Antarctic regions may be derived from photogrammetric processing of satellite images, laser/radar altimetry data, and existing cartographic data. Satellite images are a valuable source of data for obtaining a wide range of DEM products. For example, ASTER GDEM [18] from ASTER images and REMA DEM [19] from WorldView images are two representative DEMs that cover the whole Antarctica. However, there is a lack of DEMs before the launch of altimetry and 3-D mapping satellites in the 1990s. This data scarcity limits our ability to quantify surface changes in ice shelves and glaciers over long periods. Furthermore, it is difficult to carefully understand the evolution of glacier-ice shelf systems, as well as to model and predict future major calving events and even collapses.
A collection of photographs taken by United States (US) reconnaissance satellites between the 1960s and the 1970s was declassified in 1995 [20]. ARGON satellite photographs from 1963 cover the entire Antarctica and provide a broader perspective for studying the region's early conditions. These photographs also enable us to determine the surface changes that occurred in the Larsen B region during their early stages [21]. In this article, we used rigorous photogrammetric processing techniques to reconstruct a DEM of the Larsen B region from the stereo ARGON optical photograph in 1963, including its accumulation area and tributary glaciers [22], [23]. We compared the surface elevation changes between the ARGON DEM of 1963 and other DEM obtained from recent 3-D spaceborne imaging and altimetry data recorded in the 2000s to detect surface changes over 40 years before the Larsen B collapse.

A. Study Area
The Larsen B region is defined in this article as the portion of the Antarctic Peninsula between Robertson Island (65°S) and Jason Peninsula (66°S), as shown in Fig. 1 [3], [4]. The Hektoria, Green, Evans, Crane, Flask, and Leppard glaciers (from north to south) are the tributaries of the Larsen B embayment focused on in this study. The LBIS abruptly collapsed from February 2002 to March 2002, and the remnant ice shelf is known as the Scar Inlet Ice Shelf. After that, the ice mass of these tributary glaciers in the north, including Hektoria, Green, Evans, and Crane glaciers, flew directly into the Weddell Sea, while the glaciers in the remaining southern part flew onto the Scar Inlet Ice Shelf. Before the collapse in 2002, some studies were conducted using altimetry data and satellite images to determine the extent, elevation, and dynamic changes of the ice shelf. Before the collapse, there was no obvious change in ice flow velocity in the tributary glaciers and no research on the observed elevation change as early as 2001 [9], [13]. The elevation change of Larsen B tributary glaciers is tracked back to 1963 using the DEM obtained from the ARGON stereo pair in this article.

B. Research Data
The research data used in this article include satellite stereo image pairs, image mosaics, DEM products, and altimetry data (see Table I). DEM reconstruction is carried out by using The ARGON photographs taken in 1963 were initially recorded on film. They were scanned at a ground resolution of 33 m after declassification in 1995, despite the film's nominal ground sampling resolution of 140 m [23]. The scanned photographs should be able to resolve features with sizes ranging from the scanned resolution of 33 m to the Nyquist-theoryderived resolution of 59 m based on the grayscale continuity of the film [23]. The quality of the scanned photographs is relatively poor due to the imaging conditions, long-term storage of the film, and analog-to-digital scanning processes [24], [25]. As a result, additional preprocessing was required for these images. Many photographs from the three ARGON missions, 9034A, 9058A, and 9059A, were cloud-affected, especially those over West Antarctica. Missions 9034A and 9059A show noticeably lower image quality than mission 9058A. Two photographs from mission 9058A, DS09058A014MC114, and DS09058A014MC115 (see Fig. 2), were chosen to form a stereo pair for a 3-D photogrammetric reconstruction. These photographs were taken on August 29, 1963 and covered the entire study area. The photographs clearly show ice shelves, glaciers, outcrops, and mountain ridges, but they have no precise geolocation and must be georeferenced into a geodetic reference system.
The laser altimetry data and DEM products are used to provide spatial reference and analyze the elevation changes (see Table I). The WGS84 ellipsoid and Antarctic Polar Stereographic projection are used as a reference for all elevation data and products. Because no in situ data from the 1960s are available, the horizontal control was a 15-m resolution LIMA [26].
REMA DEM is a recently released DEM product that covers the entire Antarctica and is obtained from WorldView-1, 2, 3, and GeoEye-1 images from 2009 to 2017. It has an absolute accuracy of approximately 1 m and relative accuracy of decimeters [19]. ASTER GDEM and the 100-m DEM of the Antarctic Peninsula were produced using a combination of images from 2000 to 2010 and cannot provide separate terrain information prior to the 2002 collapse [18], [28]. The accuracy of the AST14DEM generated using ASTER Level-1a images and the new LP DAAC system is better than 25 m. However, there are a large number of visible blunders in relatively flat areas due to a lack of image texture and difficulties in image matching [28], [29], [30]. Thus, none of these DEM products are suitable for analyzing the elevation changes from 1963 (ARGON DEM) to 2002 (LBIS collapse) in terms of temporal coverage and data quality. Consequently, we generated a new 30-m-resolution DEM from ASTER Level 1a images, that has limited coverage and primarily includes tributary glaciers for the glacier-lowering estimates. The adopted ASTER Level 1a images included three images taken on December 22, 2001, before the 2002 collapse, that covered the majority of the glaciers; and the other three images were taken on December 7, 2002 (i.e., after the collapse), which mainly covered the accumulation area. The DEM is generated from stereo ASTER images and rational polynomial coefficient parameters using the ERDAS IMAGINE's photogrammetric processing system LPS [31]. Since no ground control points (GCPs) are used, the generated ASTER DEM must be registered with laser altimetry data, which are also used to register the ARGON DEM.
The ICESat-2 Advanced Topographical Laser Altimeter System has six laser beams arranged in three pairs, with a 90-m interval between each pair and a 3.3-km interval between beam pairs [32]. The ICESat-2 ATL06 product has a spacing of 20 m along-track and an absolute elevation accuracy of 2-4 cm [33], [34]. ICESat-2 data from March 29 to September 19, 2019, are being used to register ARGON and ASTER DEM. The study area was also covered by Pre-IceBridge Airborne Topographic Mapper (ATM) L2 laser altimetry data from 2008, with an elevation accuracy of 10 cm [35], [36], [37]. These datasets were used to test the accuracy of the reconstructed ARGON DEM in stable areas and to examine glacier elevation change.

III. PHOTOGRAMMETRIC DATA PROCESSING FOR ARGON DEM RECONSTRUCTION
The photogrammetric data processing method used to reconstruct the DEM from the two ARGON photographs is described in this section. Preprocessing was required to improve the image quality and to allow these early satellite photos to be further processed in the subsequent DEM reconstruction process proposed by Li et al. [23] for the Antarctic ice-sheet surface modelling.

A. Preprocessing
A preliminary analysis of the selected ARGON photographs revealed the presence of several data processing issues, such as speckle noise, artifacts from long-term storage, and film deformation. Furthermore, the snow and ice cover on gentle terrain caused image saturation and relatively homogeneous texture, limiting the efficiency of the standard image-matching techniques required to perform basic photogrammetric tasks, such as image orientation and 3-D reconstruction. Image enhancement was, thus, implemented during the preprocessing stage. In previous denoising tests on ARGON photographs, a mean smoothing filtering was used to denoise speckles and artifacts [38]. The adaptive histogram equalization and Wallis filtering were then used to improve the local surface texture [39], [40], [41], [42]. The film deformation was corrected during the interior orientation stage, as described in the following section.

B. Image Orientation
The image orientation stage was divided into two parts: interior orientation (IO) and exterior orientation (EO). Because only two photos were involved, these two steps were applied separately [43].
During the IO, the fiducial marks in the enhanced images were automatically recognized and used to instantiate the camera image reference system. The relationship between the measurable "image" coordinate system and the "'photo" coordinate system was established using these fiducial marks through a secondorder polynomial transformation. This IO also incorporated the focal length and lens distortion parameters from an ARGON camera calibration report. This information was also used to correct the film deformation [38].
Due to a lack of in situ data to be used as GCPs, some stable features from the ARGON stereo photographs were selected, and their coordinates were measured from the corresponding features in the available Antarctic mapping products. Peaks distributed on both sides of the glacier, as well as outcrops on the peninsula and island, were used as GCPs in the study area. These features were generally located in steep terrain, where snow is difficult to accumulate and erode over time. Despite differences in satellite sensors and illumination conditions, the same feature can be obtained by selecting the local peak from historical and modern photographs. GCPs were determined in this manner by using the LIMA mosaic with a spatial resolution of 15 m [26], and the REMA DEM with a spatial resolution of 8 m [19]. The procedure was carried out following the guidelines proposed by Ye et al. [38] and Feng et al. [44] to make sure that the GCP positions have not changed over the time span between the AR-GON images and control data. For example, outcrop peaks were selected and measured during the measurement process using orthoimages and 3-D visualization tools. We manually measured them on ARGON photographs to obtain "image" coordinates and on LIMA mosaic and REMA DEM to obtain horizontal coordinates and elevation on the geodetic reference system (see Fig. 2). Seven GCPs were evenly distributed throughout the study area and were not too close to the image borders to avoid image deformation. A rock outcrop dataset derived from Landsat 8 images [45] was used to assist in the GCP selection procedure that guarantees that the GCPs are located at outcrop peaks and no changes occurred over the ∼60 years between the ARGON and recent satellite image acquisitions (see Figs. 10 and 11). First, the rock outcrop dataset was used to determine a GCP candidate region. Then, the REMA DEM with an exaggeration factor of five was used to find the outcrop peak in the candidate region. Furthermore, we used the REMA DEM to produce a shaded relief map with the same incident azimuth and elevation angles of the sun as the Landsat 8, which is applied for true and false color interpretation of the outcrop. The Landsat 8 image was also draped on the DEM for further examination. Consequently, we successfully selected 6 GCPs on outcrop peaks. To achieve an even GCP distribution, an additional GCP, G3 (see Fig. 10), was added in a region not listed in the outcrop dataset. It is located at a sharp peak point with an ice flow velocity of less than 4.7 m/yr. The visual examination procedure described previously lead to an interpretive conclusion that the point be stable and not susceptible to snow accumulation. Thus, the point is selected as a GCP. Eight check points (CPs) were all selected on outcrops in the same way for the purpose of accuracy verification of the reconstructed orientation parameters.
To ensure the EO's stability, additional corresponding tie points (TPs) were measured on both images in the overlapping area (see Fig. 2). These TPs only required the image coordinates of the two ARGON photos' corresponding features. As a result, selecting TPs is less strict than selecting GCPs, and the number of TPs (36) is significantly greater than GCPs (6). The acquisition time between two ARGON photographs was approximately 22 s as interpolated by ephemeris file. In such a short time, the features on the two images can be considered stationary to each other. Because of accurate and clear identification of features on both images, such as peaks, ridges, ice shelves, and fast ice, the 36 TPs are evenly distributed in and around the study area. All available corresponding points (both TPs and GCPs) were used in a bundle adjustment [23], [43] that started with approximations for EO parameters and ground coordinates for TPs. The bundle adjustment allows ARGON photographs to be mapped to the surface of the 3-D ice surface.

C. Multilevel Dense Image Matching for 3-D Reconstruction
Due to the low quality and texture of the ARGON photos, standard state-of-the-art image matching techniques produced insufficient results. To address this issue, a more robust "multilevel image matching" procedure was used to combine featurebased and area-based image matching techniques hierarchically.
First, stereo-pairs of epipolar images were generated from the original scanned ARGON photographs. From the EO parameter computed after bundle adjustment, the fundamental matrix of each stereo-pair was calculated and then used to normalize each ARGON photograph [46]. The resolution of epipolar images was the identical to that of a scanned ARGON photographs (i.e., 33 m). We created a five-level image pyramid from epipolar images using Gaussian filtering. Because each level was decimated by a factor of two, the upper image (Level 4) was 1/16 the size of the original epipolar image (Level 0).
Starting from Level 4, the SIFT operator was used to detect corresponding feature points on the stereo pair [47]. This operator's robustness in finding correspondences motivated its use. In addition to the standard techniques used to examine outliers in feature-based matching with SIFT-like algorithms, these corresponding points were checked by human interpretation (e.g., [48]). Furthermore, some outcrops, mountain ridges, and glacier surfaces were manually measured and added to depict the outline of critical terrain topography. As a constraint for further image matching steps, this set of "SIFT" and "manual" corresponding points was used as seed points to generate a triangulated irregular network (TIN) in the image pyramid. After the detection of new features in both images using a Shi-Tomasi corner detector [49], feature-to-feature normalized correlation coefficient (NCC) matching was performed with the parallax and TIN constraints. Specifically, TIN was built with seed points and used to predict matches. The correspondence points on the epipolar images are distributed along the corresponding epipolar line. We set the X direction for epipolar line and the Y direction for the perpendicular line. The coordinates on the left and right images for one pair of corresponding seed points are (X 1 , Y 1 ) and (X 2 , Y 2 ). (X 1 , Y 1 , p) were used to build TIN, where p is parallax and equals X 2 minus X 1 . The coordinates (X i , Y i ) are used to interpolate parallax (p i ) under the TIN for feature points or grid points to be matched on the left image. The coordinates on the right image are then predicted as (X i +p i , Y i ). If the differences between the results of NCC and the results predicted by TIN are less than the threshold, the corresponding points are confirmed as new matches. Subpixel matching can be achieved when, for example, combined with the parabola fitting method. A surface fitting was also used to eliminate the mismatched points. Then, all seed and matched points were considered to be corresponding points. Because one pixel of the current layer represents 2×2 pixels of information on the next layer, the corresponding point must be rematched when it is transmitted to the next layer as a seed point to obtain the correct position. The parallax constraint is used during the process. Existing matching points may be lost during the rematch process due to texture information differences on each layer, especially on flat glaciers and ridges. As a result, if the distribution of seed points in the local area is poor, additional points need to be added manually.
The abovementioned feature-matching process was used at each level of the image pyramid. Finally, using the subpixel NCC matching technique, the dense matching process was performed under the constraint of TIN built by feature matches on the original size epipolar images with a grid spacing of 5 pixels [23], [38]. Although the multilayer matching method yields good results, some ridges, and glaciers with few textures do not have enough matching points. This was attributed to a variety of factors, including the fact that the shadows of the mountain peaks were elongated and affected some narrow valleys during image acquisition due to the low solar altitude angle.
The innovativeness of the abovementioned preprocessing and hierarchical matching methods is demonstrated by the effectiveness in handling the oldest satellite images and the unique glaciological conditions in the Antarctic Peninsula. First, the preprocessing method can eliminate noise and distortion caused by long-term storage and scanning of historical films. Second, the hierarchical matching method solved the challenging DEM generation problem in the unique Antarctic Peninsula environment where extremely large parallaxes differences are produced by high mountain ridges, low glacier valleys, and steep glacier margins, all within a small geographic extent. At the top layer of the image pyramid, we placed seed points at these difficult locations where regular matching methods often fail. In the following layers, new matched points can be found and propagated hierarchically to produce the full-scale DEM.

D. Quality Assessment of Image Matching and DEM Generation
By comparing distances between the matched points or features and those selected by human operators at each image pyramid layer, we evaluate the quality of the image matching results [22]. We checked for quality at each level by dividing the mapping area into grid cells. Grid cell size was determined using the scale of images. Each grid cell received one check point. From Level 0 to Level 4, 50 check points were chosen randomly throughout the study area. Furthermore, on the final 5-pixel dense-matching grid layer, 50 check points were extracted from each of the five terrain types, i.e., ice flow, peak, glacier, peninsula, and ridge (see Fig. 3). More details on this implementation strategy can be found in. The image matching evaluation result is summarized in Table II. Thus, the peaks had the lowest mean residual of 0.45±0.23 pixels. Similarly, the CPs chosen in the Jason Peninsula were generally distinct features with a lower mean residual of 0.55±0.28 pixels. A mean residual of 0.50±0.23 pixels was obtained in the ice flow area near the grounding line, where ice-flow line details and the combined line structures were used to select CPs. The glaciers had a larger mean residual of 0.61±0.31 pixels due to shadows cast from the side mountains as well as the dynamic surface. Finally, the long-ridge area on high-accumulation land was relatively flat and smooth with fewer distinct features for high-quality matching. In this

IV. RECONSTRUCTED ARGON DEM
The DEM was constructed using a total of 707929 3-D points derived from feature matching and dense grid matching. Overall,  The final matched feature and grid points were projected onto the 3-D ground coordinate system using the refined EO parameters obtained from the bundle adjustment. Although the photographs have a resolution of 33 m, which was used for image matching and 3-D positioning, the ARGON DEM was interpolated to 500 m using Kriging. This solution was motivated by the relatively low density of matched feature points in some local areas where the matching results were affected by the poor image quality and blurry texture. Fig. 4 depicts the reconstructed DEM draped with one of the ARGON photographs taken in 1963. Overall, the reconstructed terrain topography corresponds to the Antarctic surface features and imaging geometry in the study area, which includes the upper accumulation regions, glaciers, outflow areas, shadows, and relatively flat ice shelves.

A. Quality Assessment of the EO Parameters
To ensure stability to the EO computation, seven GCPs were used within the bundle adjustment, along with 36 well-distributed TPs (see Fig. 2). The numbers of GCPs and TPs were sufficient to ensure that the observations had an acceptable level of inner reliability [51]. The unit weight standard deviation was 0.52 pixels, indicating a high level of internal accuracy.
The bundle adjustment's absolute accuracy was estimated using 8 CPs (see Fig. 2). The image coordinates of the CPs were projected to the 3-D geodetic coordinate system using the estimated EO parameters. The differences between the projected ground coordinates and the positions from the LIMA mosaic and REMA DEM, obtained by manual measurement, were used to represent the errors in the EO parameters. The accuracy in X, Y,

B. DEM Registration
Both ARGON and ASTER DEMs generated in this article require registration using an accurate altimetry data set to estimate elevation changes. The horizontal and vertical shifts in the DEMs were corrected using a linear coregistration method as where X, Y are the coordinates of the DEM point and dz is the elevation difference between the DEM and the ICESat-2 points. We select a set of GCPs to estimate the coefficients a 0 , a 1 , and a 2 , each set for one DEM that are then used to compute elevation adjustments dz for the DEM points that are further    Table IV) are illustrated in shade.
used to detect changes. Feng et al. [44] proposed a systematic method for selecting GCPs using historical and modern data sets in the Antarctic environment. The GCPs are selected as points on stable landmarks, such as outcrops, peaks, blue ice, and ice rise. If these landmarks are not available, stable points with low velocities (< 10 m/yr) can also be selected. Fig. 5(a) illustrates 32 GCPs that were used to register the ASTER DEM, which consists of six sub-DEMs, with ICESat-2. The GCPs were mostly outcrops, peaks, and mountain ridges separating glaciers. Some GCPs in nonglacier areas are selected as local high locations with a velocity of less than 10 m/yr. Two Antarctic velocity maps, InSAR-based velocity map [52] and Landsat 8-based velocity map [53], were employed during the process to identify areas of low ice velocity [see Fig. 5(b)]. One set of registration coefficients are estimated for each sub-DEM. We used the elevation estimated for each sub-DEM. We used the elevation difference of the GCPs between the ASTER DEM and ICESat-2 data to show the effectiveness of registration. The overall elevation difference between the ASTER DEM from ICESat-2 data decreased from 37.2±30.9 m before registration to 0.2±2.4 m after registration.
The same registration model of (1) was also used for the ARGON-DEM. The method for GCP selection for this registration is same as that for photogrammetric bundle adjustment. A total of 9 GCPs were chosen to perform the registration of the ARGON-DEM that covers the whole region with the ICESat-2 data [see Figs. 5(c), 12, and 13]. All GCPs are selected on outcrop peaks. The overall elevation difference at the GCPs between the ARGON DEM and ICESat-2 data reduced from −2.3±6.9 m before registration to 0.0±6.1 m after registration.
The uncertainty of the ARGON DEM mainly includes two sources, matching, and registration errors. The matching error of 0.55 pixels is estimated using the checking points in the study area. Its effect on elevation estimation is then calculated as σ match for each glacier (see Table IV) based on the glacier location, ARGON mission camera, and flight information [23], [54]. Given the registration error of σ regist = 6.1 m the elevation uncertainty of the ARGON DEM is estimated as σ ARGON for each tributary glacier (see Table IV), using the following equation: Similarly, we assume that our matching error for the ASTER DEM is 0.5 pixels. Its effect on the elevation is calculated (see Table IV) m based on the same ratio above calculated for the ARGON images. Given a registration error of 2.4 m, the elevation uncertainty of the ASTER DEM is estimated for each tributary glacier (see Table IV).

A. Elevation Change in the Region and Along the Ridge at Graham Land
Many studies have confirmed that ice-shelf disintegrations cause glacier acceleration and, as a result, upstream glacier thinning. However, the changes of about 40 years prior to Larsen B's calving in 2002 had not been studied before. The 2001 ASTER DEM generated in this article covered the main tributary glacier area and was used to estimate elevation changes since 1963 (ARGON DEM). The mean difference between the two DEMs was estimated to be −6.83±91.22 m, and the R 2 of linear correlation was 0.98 (see Fig. 6). The mean difference is small, indicating that the overall differences between the ARGON DEM and ASTER DEM are insignificant. Meanwhile, the R 2 indicates the two DEMs are well consistent. Graham Land, in the northern part of the Antarctic Peninsula, is the primary mass accumulation area, supplying ice masses to the tributary glaciers of the LBIS (see Fig. 1). ATM data in Antarctica was consistent with ICESat-2 data, with differences ranging from −13 cm to 3 cm [55]. ATM data can, then, be compared with ARGON DEM after registration by using ICESat-2 data to detect elevation changes. Along Graham Land's elevation profile AA' [see Fig. 7 [56], [57], [58], snow accumulation may have increased, which could compensate for surface melting [59]. Each point in Fig. 7(a) shows the mean surface mass balance (SMB) from 1979 to 2008 using the RACMO 2.3p2 simulation results with a 5.5-km spacing [60]. Fig. 7(c) shows the average SMB within a distance of 5 km from profile AA' from 1979 to 2008. Except for a clear difference of SMB between the two sides separated by the profile AA', there is no obvious temporal change trend along the profile.

B. Elevation Change of Main Tributary Glaciers
Prior to the collapse, we used the ARGON DEM and ASTER DEM to compare the elevations of six major tributaries, the Hektoria, Green, Evans, Crane, Flask, and Leppard glaciers (from north to south), which feed the ice shelf [see Fig. 8(a)]. We calculated the post-collapse changes in the glaciers and the to the ARGON DEM from 1963. Fig. 9 depicts the elevation profiles of the six glaciers. Table V summarizes the average elevation changes of the glaciers found in this study  as well as those published in a previous study (2001)(2002)(2003)(2004)(2005)(2006) [13].
ASTER DEM and ARGON DEM have been registered by using ICESat-2 data. The ASTER DEM is of good quality and almost covers the entire glacier, from upstream to the grounding line. Despite being affected by shadows and having fewer textures in some areas, the ARGON DEM has reconstructed the terrain of the Larsen B area of the 1960s in an unprecedented manner. The ATM data (2008) [13], are only used as a reference (see Figs. 8 and 9) in this article and will not be specifically analyzed. We focus on the differences between the ARGON DEM (1963) and the ASTER DEM (2001).
Elevation differences between the ASTER DEM and ARGON DEM of these glaciers are shown in more detail in Fig. 9. The Hektoria, Green, and Evans Glaciers are located in the region far north and flow to the same subembayment. The elevation of ASTER DEM is roughly the same as that of the ARGON DEM along the Hektoria, Green, and Evans Glaciers [see Fig. 9 Overall, the elevation change of all six tributary glaciers are all within the elevation uncertainties (1σ) and, thus, insignificant. The estimated elevation change rates are also insignificant in comparison to those after the 2002 collapse (see Table V). As shown in Fig. 9, the Hektoria, Green, Evans, and Crane Glaciers did not change significantly before the 2002 collapse, but declined rapidly afterward, showing an elevation decrease rate range from less than 1 m/yr (precollapse) to nearly 20 m/yr (postcollapse).
We show that there are no elevation changes for the two southernmost glaciers, Flask and Leppard Glaciers. After the collapse, both glaciers were still supported by the ice-shelf remnant. Therefore, considering the accuracy of both DEMs, Fig. 13. Illustration of multisource data interpretation used in the rigorous GCP selection procedure, from Point GR1 (a) to Point GR9 (i). Within the panel of each point, (i) is the Landsat 8 color image; (ii) shows the outcrop (brown polygon) areas from a dataset produced using Landsat 8 multispectral images [45]; the GCPs are located on peaks of outcrops of the dataset (except G3); (iii) and (iv) are Landsat 8 and ARGON images overlaid on the REMA DEM, respectively; the viewing direction is defined in (i).
there is no significant elevation change on the six glaciers from 1963 to 2001.
In Section V-A, we observed that there is no obvious elevation change in the accumulation area (around line AA') and no change trend in the snow input estimated from the RACMO 2.3p2 model. Moreover, there are no velocity accelerations from 1996 to 2000 by analyzing the velocity maps extracted from ERS-1/2 images [9]. This indicates that the input and output mass of these glaciers were relatively stable and balanced between 1963 and 2001. It could provide evidence that the ice shelf supported the inland glaciers over nearly four decades before the collapse. The elevation of these glaciers did not drop significantly until the collapse occurred. The results demonstrated that the hierarchical image matching method can be modified and applied to reconstruct a historical Antarctic DEM using satellite images acquired ∼60 years ago through an innovative and rigorous ground control point selection procedure that guarantees no changes occurred at these points over the period. The average accuracy of the ARGON-derived DEM and ASTER-derived DEM are ±42.87 m and ±17.70 m, respectively. Our results reveal that the six major tributary glaciers in the embayment, Hektoria, Green, Evans, Crane, Flask, and Leppard Glaciers, had insignificant elevation changes before the collapse. In addition, there are no significant elevation changes in the Graham Land ice mass supply area both before and after the collapse. These findings are consistent with the trend that there was no increase in mass loss and no acceleration of ice flow in the upstream basins between 1996 and 2000 [9]. The combined analysis of the newly produced ARGON DEM in 1963 and ASTER DEM in 2001 suggests that the Larsen B region tributary glaciers were in relative equilibrium before the 2002 collapse. These DEMs fill a decades-long elevation gap in the survey She is currently a Lecturer with the Department of Surveying and Geo-Informatics, Tongji University, China. Her main research interests and education responsibilities include photogrammetry, remote sensing, and applying techniques and datasets to monitor and understand mass and energy changes of glaciers and ice sheets. Her research direction is surveying and mapping, specialized in photogrammetry with multisource remote sensing data. Her research interests include the reconstruction of DEM from multisource satellite images and instability of ice shelf.
Lu An received the Ph.D. degree in earth system science from the University of California Irvine, Irvine, CA, USA, in 2017.
She is currently an Associate Professor with the College of Surveying and Geo-informatics, Tongji University, Shanghai, China. She has worked to use high-resolution airborne gravity data in combination with other datasets to determine the bathymetry of fjords and the bed topography of glaciers in Greenland and Antarctica. These datasets are critical to understanding the role of ice-ocean interaction in controlling the evolution of glaciers and ice sheets under global warming. Her research interests focus on better understanding and explaining ongoing changes in the cryosphere and reducing uncertainties in the contributions of ice sheets to sea level rise using various satellite data. Marco Scaioni received the M.Sc. degree in civil engineering and the Ph.D.