Displacement estimation from image matching in the spatial domain

The foundation for displacement estimation from repeat remote sensing is pattern matching (Scambos and others, Reference Scambos, Dutkiewicz, Wilson and Bindschadler1992; Kääb and Vollmer, Reference Kääb and Vollmer2000; Strozzi and others, Reference Strozzi, Luckman, Murray, Wegmüller and Werner2002; Debella-Gilo and Kääb, Reference Debella-Gilo and Kääb2011). It is based on finding a similar template from an image in another larger subsection of an image at a different time instance (t _q). Often the similarity measure (ρ) is implemented through normalized cross-correlation as,

(1)$$\rho_{\,pq} = {\sum_{k\comma l}{\lpar {\bf I}_{\,p}\lpar i + k\comma\; j + l\rpar -\bar{{\bf I}}_{\,p}\rpar \lpar {\bf I}_{q}\lpar k\comma\; l\rpar -\bar{{\bf I}}_{q}\rpar }\over \sqrt{\sum_{k\comma l}{\lpar {\bf I}_{\,p}\lpar i + k\comma\; j + l\rpar -\bar{{\bf I}}_{\,p}\rpar ^{2}} \sum_{k\comma l}{\lpar {\bf I}_{q}\lpar k\comma\; l\rpar -\bar{{\bf I}}_{q}\rpar ^{2}}}}.$$

Here, ρ_pq denotes a score surface at discrete step sizes (k, l). It is calculated by the similarity of two templates (I_p and I_q) at different time instances (t _p and t _q). The upper part of the fraction of equation (1) shifts the intensities in the templates by their mean intensity (denoted by the overbar) to transform the signals to zero-centered variations, while the lower part of the fraction normalizes the contrast of the templates (or intensity spread). This results in a similarity score that can range between −1 and +1.

Matching the velocity over an icefall, or other similarly fast changing surfaces, is much more challenging than on a slow and homogenous flowing glacier tongue. Surface features in an icefall can change quickly over time, thus visual distinctive features, which can be tracked over time, are quickly lost. Furthermore, the flow is first converging and then diverging and thus exhibits considerable shear. Consequently, bulk velocity estimation (of homogenous, simple translation) does not (or not well) work, and in order to fulfill this constrain smaller template sizes are needed, which then lowers the signal-to-noise level (Debella-Gilo and Kääb, Reference Debella-Gilo and Kääb2012). Therefore, in this specific case, we apply pattern matching in the time domain (Meinhart and others, Reference Meinhart, Wereley and Santiago2000), known as ensemble matching (Delnoij and others, Reference Delnoij, Westerweel, Deen, Kuipers and van Swaaij1999), see also the schematic in Figure 2. This means the matching is done with much smaller windows than necessary for common single image pair matching, but instead information is gathered from multiple image pairs of a whole collection. Because some of the imagery will have cloud cover, not all image pairs can be used from the stack, and therefore a selection is made and this we call a collection. With ensemble matching, the same or higher signal-to-noise can be achieved, as a similar total number of image samples (pixels) is taken as is done with single image pair matching, i.e. replacing spatial redundancy by temporal redundancy. When the flow does not change over the analyzed collection period, a clear peak will emerge. Laboratory experiments of sparsely illuminated seeds within a fluid or airmass have shown this spatial resolution can be reduced to a single pixel resolution of independent flow estimates, when an image stack size of several hundreds or thousands of images is included (Westerweel and others, Reference Westerweel, Geelhoed and Lindken2004). For a highly fractured crevasse field, a smaller sample might suffice.

Fig. 2. Schematic illustration of the difference between regular single image pair matching and time averaged ensemble matching. Left: single image pair matching requires the application of large window sizes to achieve smooth and clearly defined correlation maxima. Middle: The application of small window sizes and single image pair matching leads to noisy correlation surfaces with small signal-to-noise ratios. Right: Within ensemble matching, the individual noisy correlation surfaces are combined (here: summed) over time to achieve clear correlation peak with high signal-to-noise ratio despite small window sizes.

Sub-pixel localization

Several methods exist to estimate the localization of the correlation peak, which is inherently at pixel level, at sub-pixel resolution. A main difference between the method's implementation is the amount of information needed from the correlation scores surrounding the maximum. Either the directly neighboring pixels are used that are next to the highest score, independently estimating the peak along each axis. A second feature to consider is the shape of the peak, which can be modeled by a parabolic (Argyriou and Vlachos, Reference Argyriou and Vlachos2005), Gaussian (Willert and Gharib, Reference Willert and Gharib1991) or triangular (Olson and Coombs, Reference Olson and Coombs1991) function. Other approaches use all pixels surrounding the highest correlation score, which is the case for center-of-mass, a 2D Gaussian (Nobach and Honkanen, Reference Nobach and Honkanen2005), or a spline (Fahnestock and others, Reference Fahnestock2016) function. The former methods are most popular due to their low computational power, while more computationally expensive methods like upsampling of the imagery and/or its correlation surface is also possible and can achieve precisions in the order of 1/20 of a pixel (Debella-Gilo and Kääb, Reference Debella-Gilo and Kääb2011). In this study, we investigate all the above mentioned methods against upsampling, as systematic errors are present in the different localization methods. It is therefore of importance to assess these different methods for their preferences to integer value displacements, known as peak-locking (or pixel-locking).

Sentinel-2 high repeat satellite imagery

The EU/ESA Copernicus Sentinel-2 mission consists of a constellation of two medium resolution satellites recording visible to shortwave infrared imagery with a spatial resolution down to 10 meter. The first satellite of this mission, Sentinel-2A, was launched in June 2015 and has a same-orbit repeat cycle of 10 days. With the launch of Sentinel-2B in March 2017 into the same orbital plane, this revisit decreased to 5 days from the same orbit. The high dynamic range (12 bit raw) of the Sentinel-2 images makes them ideal for glacier velocity extraction as also areas with little surface contrast show often patterns useful for offset tracking. For correct displacement estimates, imagery from the same orbital track needs to be used (Kääb and others, Reference Kääb2016), as orthorectification errors due to erroneous or outdated elevation models distort the images. Methods are available to correct for these artifacts and make the use of imagery from alongside orbits possible (Altena and Kääb, Reference Altena and Kääb2017). However, in this study, such methods are not the focus and not needed, as the consistent same-orbit monitoring by Sentinel-2 provides enough imagery suited to demonstrate the method and derive a useful velocity field over the Khumbu icefall.

Location-aware sport wearables

For independent comparison to our velocity data, we use trajectory data collected by high-altitude mountaineers. These climbers move through the Khumbu icefall while equipped with location-aware wearables. These handheld GPSs have sensors that are able to record position, altitude and heading during outdoor activity. Afterwards the data can be uploaded to online sports diaries, where these can be shared either within the respective communities or archived into a private collection. Accumulation of such spatio-temporal datasets can be used to extract geographic information, mostly done through trajectory clustering, patterns or prediction (Giannotti and Pedreschi, Reference Giannotti and Pedreschi2008).

These types of positioning data can be used in different ways, depending on spatial and temporal selection criteria. Definitions are therefore given about such groupings for this data-type. For instance, a ‘trajectory’ is a collection of positioning points connected to each other by an increasing and ordered timestamp. A ‘trail’ is a collection of positioning points on a line without a time instance. If several trajectories are close together in space and time they are called a ‘flock’. Then, when a flock is in close proximity to a trail, they are called a ‘convoy’.

During a climbing ascent through the icefall, the data collection is thus as a convoy. However, in our case, the multi-temporal sequences can be seen as being within an environment that changes over time (Wiratma and others, Reference Wiratma, van Kreveld and Löffler2017). In our situation, we assume the flock follows a trail, but the trail has a temporal character; as the glacier moves with the trail on top. Consequently, glacier velocity can be extracted over time from different flocks.

Cloud-free scene classification

Cloud cover hampers the ability to extract within-season imagery for many optical Earth observation systems like Landsat (Ju and Roy, Reference Ju and Roy2008). The optical instruments of Sentinel-2 are also hampered by the cloud coverage, but this constellation has an increased revisit rate of 5 days, increasing the possibility to get many pairs of images with a short acquisition interval even if a number of scenes are affected by cloud cover. Due to the resulting large data volumes, the collection of useful imagery requires automatic classification of cloud cover. Global cloud cover products in the Sentinel-2 meta-data can be used, but these would exclude partly clouded or highly clouded scenes, even if the very location of interest (here the icefall) is actually cloudfree. The detection of clouds in high mountain terrain is challenging due to snow cover and topographic shadowing. Therefore, we use all the images and stack statistics (median) at the pixel level to generate a scene composite (or stack-median, $\tilde {{\bf I}}$) (Winsvold and others, Reference Winsvold, Kääb and Nuth2016). The stack-median is then converted by a 2D Fourier transform, and its absolute values are taken and transformed by their natural logarithm ($\ln {\Vert {\cal F}\lpar \tilde {{\bf I}}\rpar \Vert }$). These operations result in a distinct pattern description, particularly present in natural images (Torralba and Oliva, Reference Torralba and Oliva2003), and gives a mean to compare against other transformed images. Thus, every image is transformed by these operations, and the resulting image description is matched against the stack-median. The similarity metric calculated for this comparison is mutual information (Viola and Wells III, Reference Viola and Wells III1997). These mutual information scores are ultimately used to classify imagery into cloud-free and cloud-covered. This class separation is achieved through a two-class k-means procedure, resulting in an binary cloud-cover attribute to each image that is used to build the collection of images suitable for ensemble matching.

Icefall displacement from Sentinel-2

Currently, the absolute localization of Sentinel-2 is done through its sensor readings, hence its variability is in the order of 10 meter (Gascon and others, Reference Gascon2017). Each image is therefore in our study co-registered by simple translation to the stack-median, using conventional image matching over stable terrain (Leprince and others, Reference Leprince, Barbot, Ayoub and Avouac2007), and each individual image is interpolated towards this base image. The stack-median is thus also affected by the same co-registration noise. However, these varying offsets do not seem to cause much blurring in the stack-median. This might however be the case for image collections with worse registration accuracy than Sentinel-2, or when significantly less images are stacked.

The individual images are then ordered into pairs from the same relative orbital track (in this case R076) and matched if they are separated by 10 days. This time-span seems to be a good compromise between significant spatial movement in relation to degradation of traceable appearance as ice flows through the icefall, where it melts and fractures. Hence, imagery of the same Sentinel-2 orbit and satellite is matched with each other. In this way also data from January 2016 onwards can be used, the time when only Sentinel-2A was in orbit. The last acquisitions used in this study are from October 2019. Given the time-span constrain of 10 days, a total of 53 pairs are matched against each other, with a template (I_p) size of three pixels, while the search space is five pixels wide.

Icefall displacement from GPS trajectories

The mountain climber trajectories of an ascending (3rd of May) and a descending track (18th of May) were used for comparison to our surface velocity estimate. The tracks come from a Norwegian climbing team that used Garmin inReach satellite positioning devices as they climbed through Khumbu icefall in 2019, toward Camp II and Mount Everest. One of the devices has been set to record every minute, while another member recorded at a 10 min interval (· and + in Fig. 1, respectively). For a large extent, the path goes straight against the flow direction of the glacier. In addition, the positioning data are very noisy, most presumably due to multi-path signals from the surrounding peaks and steep ice walls (i.e. the large spread in elevation in the subplots of Fig. 1). There is only a small part of the route where the trail is perpendicular to the flow, and it is this region, which is here used as validation against the ensemble matching estimate. The trajectories were first pre-processed through a Hampel filter (Lee and Krumm, Reference Lee and Krumm2011) and then a cubic smoothing spline was interpolated through these points. A smoothing parameter of 5 · 10⁻⁵ was chosen and checked for consistency with the 10 min-interval dataset. At the velocity sampling distance of 30 m, the closeness (Wiratma and others, Reference Wiratma, van Kreveld and Löffler2017) of the ascending trajectory to the descending trajectory in the direction of glacier flow was then measured to extract a glacier velocity estimate.