ICE FLOW VELOCITY MAPPING IN GREENLAND USING HISTORICAL IMAGES FROM 1960s TO 1980s: SCHEME DESIGN

Ice flow velocity over long time series in Greenland plays a vital role in estimating and predicting the mass balance of the Greenland Ice Sheet and its contribution to global sea level rise. However, there are few Greenland ice flow velocity products with large spatial coverage available showing the Greenland ice flow velocity pattern before the 1990s. We proposed three methods, including parallax decomposition, grid-based NCC image matching, feature and gird-based image matching with constraints for estimation of surface velocity in East Antarctica based on ARGON KH-5 and LANDSAT imagery, and a systematic compilation method for the ice surface velocity in East Antarctica from the 1960s to 1980s. Based on the above methods and cartographic processes, this study designed a framework for the mapping of the historical ice flow velocity of the Greenland ice sheet. Currently, the early LANDSAT images covering several glaciers in North Greenland and Northwest Greenland have been processed and applied to velocity mapping using the cartographic process proposed in this study, and some preliminary results have been obtained.


INTRODUCTION
The mass loss of the Antarctic Ice Sheet (AIS) and the Greenland Ice Sheet (GrIS) is the largest contributor to global sea level rise in recent decades (Rignot et al., 2011). The GrIS is second only to the Antarctic ice sheet in area, but it is more fragile given its recent rapid melting rate, level of mass loss, and increasing instability. Compared with the AIS, the GrIS is of lower latitude, with its southernmost end being near 59°N. The combined effects of ocean and atmospheric warming have caused an increase in ice flow and surface melt runoff, accounting for 60% of the recent contribution to the global sea level rise (10.6 mm) when compared with the AIS (7.2 mm) since 1990s (Slater et al., 2020). These phenomena reflect the global significance of the GrIS in terms of environmental change. Climatologically, it is a large regional cooling centre which plays an important role in indicating global climate change (Weidick et al., 1995). Greenhouse gases are stored in the GrIS and melting glaciers can contribute to greenhouse gas emissions (Lamarche-Gagnon et al., 2019). Hydrologically, it holds large quantities of fresh water which is an essential part of the hydrological cycle. If fully melted, it would be equivalent to a rise in global sea level of about 6.5 m (Slater et al., 2021). From the point of view of glacial hazards, it releases large amounts of icebergs into the fjords and surrounding oceans every year (Weidick et al., 1995). Therefore, the accurate measurement of Greenland glacier velocity is of great importance to studying the mass balance of GrIS and predicting future sea level rise. A long-time series of ice flow velocity maps is necessary to fully understand the evolution mechanism of the GrIS.
Currently, the published methods and ice flow velocity maps covering the GrIS are mainly as follows. Rosenau et al. (2015) combined the two feature tracking methods of fast normalized cross correlation (NCC) and least-square matching (LSM) with * Corresponding author an affine geometric model to track the features of Landsat 1-8 images from 1972 to 2015, providing a monitoring system for area-wide flow-velocity fields of the outlet glaciers of the GrIS. Fahnestock et al. (2016) used a extension of normalized crosscorrelation as applied in IMCORR to track the displacement of large-area Landsat 8 images. The produc was provided by the NASA MEaSUREs Global Land Ice Velocity Extraction (GoLIVE) project . The data set provides ice flow velocity maps with a resolution of 300 m from 2013 to 2017. The NASA MEaSUREs ITS_LIVE project (Gardner et al., 2019) provided velocity data generated using auto-RIFT (Gardner et al., 2018). These datasets used Landsat 4,5,7 and 8 imagery to generate Greenland-wide ice flow velocity products for 1985 to 2020. The resolution of the annual velocity maps is 240 m, and the resolution of the complete time average velocity map is 120 m. Mouginot et al. (2017) proposed a synergistic approach to process Landsat-8 (optical), Sentinel-1, and RADARSAT-2 (SAR) data and automatically calibrated, mosaiced, and integrated these data sets into seamless and icesheet-wide products. Mouginot et al. (2019) derive surface ice flow velocity from 1972 to 2018 using data from 16 satellite sensors deployed by 6 different space agencies (Joughin et al., 1998;Michel and Rignot, 1999;Mouginot et al., 2012), with a resolution of 150 m. Joughin et al. (2010) mapped ice flow velocity map over much of the GrIS between 2000 and 2005, using RADARSAT. Joughin et al. (2018) processed a large number of SAR and Landsat 8 images collected from 1995 to 2015 and generated a nearly complete 250-meter resolution Greenland ice flow velocity map, and they totally released 8 different products using SAR and optical images, covering the period from 1985 to 2021. These products are currently available at the US Snow and Ice Data Centre (NSIDC) and are part of the NASA Making Earth System Data Records for Use in Research Environments (MEaSUREs) program. Overall, products utilizing Landsat7-8 and SAR imagery account for most of the released products due to the high resolution and high precision positioning. Current studies on large-extent GrIS ice flow velocity have focused on the post-1985 period, and pre-1985 glacial flow velocity studies are mainly concentrated on some typical areas, such as glaciers in the north and northeast.
Satellite remote sensing images covering the GrIS date back to declassified intelligence satellite photography (DISP) imagery in the 1960s (Zhou et al., 2002) and Landsat MSS and TM images in the 1970s and 1980s. These images have not yet been fully used to systematically map the ice flow velocity at the ice sheet scale, mainly due to the low quality and complexity of image matching. Li et al. (2017a and2017b) and Ye et al. (2017) presented three methods for reconstructing surface velocity field in East Antarctica from the 1960s to the 1980s based on ARGON KH-5 and LANDSAT imagery. Li et al. (2018) presents a set of systematic maps compiling methods of developing ice surface velocity products for the entire East Antarctic area. Based on the previous studies, this paper studies and explores the feasibility for mapping ice flow velocity in Greenland using the above ice flow velocity measurement methods and velocity product preparation scheme. We present a preliminary scheme design of the historical ice flow velocity map in Greenland.

Map projections
Selecting the appropriate projection is the foundation of mapping. Projection details vary by region. Landsat imagery uses the Universal Transverse Mercator (UTM) map projection and Polar Stereographic for Antarctica (https://doi.org/10.5066/F7H994GQ). This means that images of Greenland are not uniformly projected in the same projection as those for the Antarctic region. Therefore, there is an additional step of coordinate in reference systems.
In this paper, we chose the NSIDC Sea Ice Polar Stereographic projection, which is based on the WGS 1984 ellipsoid (Northern Hemisphere: EPSG 3413), as described in Table 1. This projection with 70°N as the standard latitude mostly used in polar research, especially small-scale mapping of satellite imagery and mosaic of Greenland because of the minimal distortion around the margin of Greenland (Snyder, 1982). In addition, many Greenland data products use this coordinate system. NSIDC encourages all new products to use EPSG codes 3413 (https://nsidc.org/data/polar-stereo).
Projection conversion is involved due to applications of different image and map reference systems. We provide two solutions. One is to complete the matching in the UTM coordinate system where the image is located, and then reproject the matching results to the map coordinate system. The other is to reproject the image to the mapping coordinate system before matching and mapping. After verification, we found that the features of glaciers in the images would be distorted due to the resampling after the image was reprojected (Figure 1). In order to reduce such distortion, we decided to pre-process and match the images under their respective UTM projection zones before transforming the matching results into the WGS 84 / NSIDC Sea Ice Polar Stereographic North coordinate system.  The same features are shown in the same boxes. Although the overall resemblance is good, some pixels have changed.

Mapping scale
The glaciers in Antarctica and Greenland belong to continental glaciers, but the environments of these glaciers are quite different. The AIS is much larger than the GrIS, and most of the major glaciers are also large in size than Greenland glaciers. In the Antarctic, for large glaciers, one Landsat image can only cover part of the glacier (Figure 2b). To cover it completely, multiple images are required. For Greenland glaciers, one image can cover multiple glaciers (Figure 2a). This leads to the complexity in feature matching with a significant velocity variation in one image pair and different map compiling techniques. Here we mention two mapping scales: image scale and glacier scale. The image scale is measured in terms of the coverage of a pair of images. The glacier scale is measured in terms of the coverage of individual glaciers.
In the Antarctic region, the velocity variation expressed by a pair of images is only limited to that in a part of the glacier, that is, the change fluctuation is relatively small. For the Greenland region, the velocity change information displayed by a pair of images is significantly higher, that is, high and low velocities are scattered in one scene. This will increase the possibility of mismatching.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France In addition, due to the large area of Antarctic glaciers and the wide drainage area, the surface textures formed on the Antarctic glaciers are larger, while most of the glacier outlets in Greenland are in the fjords, with faster flow velocity and narrower passes, resulting in the smaller image texture. Limited by the quality of the early images and the constant resolution, the surface of the Greenland glaciers in the early images looks more "smoother". Therefore, in the subsequent mapping process, it is necessary to handle smaller areas within one image pair to match features of similar velocity in a "zoom in" region. To address this issue, we reduced the mapping scale from image level to glacier level. The specific procedure will be discussed in Section 3, and the comparative results of the experiments will be shown in Section 4.

FRAMEWORK DESIGN OF ICE FLOW VELOCITY MAPPING IN GREENLAND
In response to the questions raised in Section 2, we have adjusted the mapping process for the AIS previously proposed by the research group (Li et al., 2018). The following are the preliminary design and implementation steps.
Considering the topography of GrIS, the images covering Greenland and the existing velocity products, we conducted a pilot study first before a full implementation was performed. In terms of image time selection, Landsat1-5 MSS images from 1972 to 1984 are the main ones, followed by two DISP mosaics in 1962 and 1963 (Zhou et al., 2002), and finally Landsat4-5 TM images from 1985 to 1989. In terms of region selection, we process the coastal region first and then inland. We also consider large glaciers over small ones.

Image selection
When extracting ice flow velocity information from remote sensing images, it is ideal to use cloudless images and filter the images with the cloud cover option when downloading the images. However, due to the particularity of Greenland's environmental conditions, it is impossible to achieve this in many regions. Therefore, in most cases, this paper selects the best image by visually checking all candidate images. In some areas where clouds continue to appear, the measurement results obtained from the cloudless area of multiple images are compiled to obtain the best data for the area. If the data gap is small, the ice flow velocity is interpolated. In contrast, if the gap is large, it will be left blank in the final data product to ensure the high quality of the product.
For the selection of image pairs, it is important to consider not only cloudiness but also the illumination conditions and the time interval between images. As many glacial outlets in Greenland are in fjords, shadows due to differences in elevation of the terrain are unavoidable. In order to reduce mismatching caused by grey scale variations, two images with similar sun azimuth angles are selected for image matching wherever possible in the optional images. In addition, as some glaciers in Greenland have faster flow rates and shorter glacier lengths than those in Antarctica, images with intervals of several months to a year are appropriate to avoid distortion or disappearance of features.

Image pre-processing
Image pre-processing includes image enhancement and orthorectification.
In this paper, the adaptive histogram equalization method is used to enhance the images. This method optimizes the local image contrast and divides the image into regional grids for local histogram equalization. Adaptive histogram equalization with restricted contrast limits the enhancement range of local contrast on the basis of adaptive histogram equalization, which can effectively deal with noise and excessive local contrast (Muniyappan et al., 2013).
Although the Landsat section images (L1T product) were orthorectified using Ground Control Points (GCPs) and Digital Elevation Models (DEMs) from the Global Land Survey (GLS) 2000 datasets, these images still have large geometric deformation, and the positioning error in some areas can reach more than 1 km (Rosenau et al., 2015). Compared with a true orthophoto, this Landsat product still has large errors. Therefore, in order to improve the accuracy of the reconstructed velocity field, orthorectification is needed.
We used distinct feature points in bare rock areas as basic control points and Toutin's model (Toutin, 2004) in commercial software PCI to orthorectify remote sensing images (Li et al., 2018). We adopted the horizontal and vertical references of the Image Mosaic v1 and DEM from GeoEye and WorldView Imagery v1, respectively, of the MEaSUREs Greenland Ice Mapping Project (GIMP) 2000 (Howat et al., 2014).

Image matching
In this paper, we used the methods of image matching proposed by Li et al. (2017a and2017b). The hierarchical matching of feature points and grid points with additional triangulation constraints is introduced to improve the accuracy and density of matching.
Before matching, we use a mask to separate glacier regions from the background in an image. The matching parameters to be set include mainly the size of the reference window and the search radius. The size of the reference window determines how large the feature neighbourhood is to calculate the correlation. With the increase of the window, more information is contained. We use 16 or 32 pixel reference windows. The appropriate search radius is determined by the following formula (Liu et al., 2012): where r = minimum search radius vmax = the maximum ice flow velocity within the image coverage. This value is the predicted value, which can be obtained according to the existing ice flow velocity products or manual measurement points ∆t = timespan of image pair p = spatial resolution of image.

Gross error elimination
The strategy of this paper is to adopt a strict false matching elimination standard to ensure the accuracy of the matching points. The gross errors are eliminated based on the correlation coefficients of matching points, the ice flow direction, and the existing ice flow velocity products (Li et al., 2017).

Interpolation
The results of feature matching are discrete points, and the real glacier velocity field is a continuous surface, so the last step of the ice flow velocity extraction process is to use all the feature matching results to fit the ice flow velocity map through an appropriate interpolation method. In this paper, the natural neighbourhood method is selected for interpolation, and the resolution of the interpolated image is determined according to the density of matching points. As it is still in the production stage, the results are converted to the Polar Stereographic North coordinate system for interpolation. As the mapping products vary considerably from region to region, the resolution of the final product needs to be further determined based on an analysis of all regions.

EXPERIMENTS
In this paper, two different types of glacier areas are selected as examples. The locations of the glaciers are shown in Figure 3. We use Greenland basin boundary provided by Mouginot and Rignot (2019). They divide Greenland, including its peripheral glaciers and ice caps, into 260 basins grouped in seven regions: southwest (SW), central west (CW), northwest (NW), north (NO), northeast (NE), central east (CE), and southeast (SE).

Glaciers in North Greenland
Among the seven regions in Greenland, the mapping situation in the north is the closest to the East Antarctic. Glaciers account for a large proportion of the image, and the image features are clear and obvious. Figure 4 shows the matching results of the two glaciers in the north, and Figure 5 shows the final velocity map.  (1975.04-1976.04 and 1983.04-1984.05). The black boxes are the ranges of matched images. The orange and pink boxes indicate the range of the two matching periods respectively. The background is from Landsat Image Mosaic of Greenland (Chen et al., 2020).  images (1975.04-1976.04 and 1983.04 -1984.05). The basin boundary is from Mouginot and Rignot (2019). The background is from Landsat Image Mosaic of Greenland (Chen et al., 2020).

Glaciers in Northwest Greenland
The Tracy and Heilprin glaciers in Northwest Greenland are selected for comparison with glaciers in the north Greenland. Figure 6 shows our matching results at image scale and glacier scale. The number and quality of the matched points at glacier scale are better than those at image scale. Figure 7 shows the preliminary mapping results at the glacial scale.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France  . Preliminary velocity mapping results for the Tracy and Heilprin glaciers in Northwest Greenland using early Landsat images (1973.03-1973.09). The basin boundary is from Mouginot and Rignot (2019). The background is a Landsat Image Mosaic of Greenland (Chen et al., 2020).

CONCLUSIONS AND FUTURE WORK
By analysing the differences in velocity mapping between the Antarctic Ice Sheet and the Greenland Ice Sheet, and according to the compilation method of velocity mapping in East Antarctica summarized by the research group, this paper puts forward a design scheme for ice flow velocity mapping in Greenland from the 1960s to 1980s based on early DISP and Landsat images. Our goal is to produce an ice flow velocity map of Greenland from the 1960s to 1980s. We are now in the initial stage of product production. According to our experiments and exploration, we have completed the pilot velocity maps of some glaciers.
In the future, due to the terrain and image quality problems in the inland of the Greenland Ice Sheet, it is necessary to conduct further experiments to explore and improve the mapping scheme. For the accuracy evaluation of subsequent flow rate products, the error caused by projection conversion should also be considered. Further quality evaluation and scientific analysis will be done after the product is generated.