Correlation dispersion as a measure to better estimate uncertainty in remotely sensed glacier displacements

In recent years a vast amount of glacier surface velocity data from satellite imagery has emerged based on correlation between repeat images. Thereby, much emphasis has been put on the fast processing of large data volumes and products with complete spatial coverage. The metadata of such measurements are often highly simplified when the measurement precision is lumped into a single number for the whole dataset, although the error budget of image matching is in reality neither isotropic nor constant over the whole velocity field. The spread of the correlation peak of individual image offset measurements is dependent on the image structure and the non-uniform flow of the ice and is used here to extract a proxy for measurement uncertainty. A quantification of estimation error or dispersion for each individual velocity measurement can be important for the inversion of, for instance, rheology, ice thickness and/or bedrock friction. Errors in the velocity data can propagate into derived results in a complex and exaggerating way, making the outcomes very sensitive to velocity noise and outliers. Here, we present a computationally fast method to estimate the matching precision of individual displacement measurements from repeat imaging data, focusing on satellite data. The approach is based upon Gaussian fitting directly on the correlation peak and is formulated as a linear least-squares estimation, making its implementation into current pipelines straightforward. The methodology is demonstrated for Sermeq Kujalleq (Jakobshavn Isbræ), Greenland, a glacier with regions of strong shear flow and with clearly oriented crevasses, and Malaspina Glacier, Alaska. Directionality within an image seems to be the dominant factor influencing the correlation dispersion. In our cases these are crevasses and moraine bands, while a relation to differential flow, such as shear, is less pronounced on the correlation spread.

1 Introduction

The increased global availability of satellites images has created unprecedented archives of velocity products over glaciers, ice caps (Fahnestock et al., 2016; Millan et al., 2019; Friedl et al., 2021) and ice sheets (Rosenau et al., 2015; Joughin et al., 2018). These velocity fields have a large potential to enhance our understanding of ice mechanics and glacier dynamics in space and time. Current efforts are mostly focused on the automatic construction of large-scale time series (Gardner et al., 2018; Altena et al., 2019; Derkacheva et al., 2020) or the detection of special speed variations, such as seasonal fluctuations or surge dynamics, from a patchwork of velocity products (Greene et al., 2020; Riel et al., 2021). Advances in time series construction of glacier velocities will likely mature rapidly in the next few years with the new and increasing availability of suitable data. One promising direction of development is to include the measurement precision into the estimation procedure for glacier velocity variations, through either Bayesian inferences (Brinkerhoff and O'Neel, 2017) or generalized least squares (Altena and Kääb, 2017; Riel et al., 2021). Though, for such approaches estimation of the dispersion of individual image correlations is needed. Dispersion in this context is the magnitude of fluctuation or the expected variability in the velocity estimate (i.e., variance σ). Typically, a constant variance is set for the whole dataset (known as homoscedasticity), as well as an absence of correlation (ρ) between observations of different velocity components (Leprince et al., 2007). The dispersion estimation is then based upon sampling statistics, using a region of bare and stable ground if available, to extract a group variance along each axis (Herman et al., 2011; Heid and Kääb, 2012). However such bare ground might be a correct representation neither for glacier surfaces nor for their correlation dispersion estimate, as the image content and in particular the characteristics of image contrast to be matched are typically different between off- and on-glacier area and vary in addition across the glacier surface.

In our opinion the assumption of constant variance (homoscedasticity) does not hold, as displacement extraction is based upon pattern matching of small subsets of imagery, where the image content influences the displacement precision. Pattern matching is based upon similarity metrics between the matched image across its extent. Such an image subset can have texture with a strong directionality, such as crevasses, or the texture in an image subset is distorted due to skewed flow, such as shear (Debella-Gilo and Kääb, 2012). Both effects are common on glaciers but vary across the scene, thus variation in dispersion might occur across a scene as well. Within image matching, similarity between imagery is computed for a multitude of locations within a neighborhood, resulting in a surface of correlation scores for each potential displacement location, and the maximum peak of this surface is typically detected to indicate the most likely image offset. Since neighboring displacement locations have a similar appearance and partial overlap, the similarity score captures smearing in the form of an elongated spread of the correlation peak; i.e., such a peak is not a sharp spike but rather a smooth top or dome. For distinct directionality in the matched pattern location, the correlation surface gets elongated in the prevailing direction, and such effects can thus be used to extract a better formulation of dispersion for that specific matching location and time interval.

The issue of homoscedasticity can also be approached from the perspective of optical flow. Pattern matching and optical flow can be seen as interchangeable techniques, as they are mathematically similar (Lemmens, 1988). If image gradients are present in several directions, the span of the matrix is sufficiently large. However, when there is a predominant direction in the image gradients, the matrix becomes rank deficient and the optical flow estimation becomes ill-posed (also known as the aperture problem). Hence, treating the displacement axis independent does not hold, nor is a fixed precision term sufficient.

In this contribution, we demonstrate a fast estimation approach for dispersion characteristics for individual displacement estimates from image matching. These dispersion characteristics are then used to explore the connection between the correlation spread and the processes of shear flow and crevasse orientation. This gives a better understanding of the image regions where displacement estimates need to be interpreted with caution. Furthermore, our method enables better quantification of error propagation into the remote sensing and derived model products, which can improve inferences about, for instance, strain rates, glacier depth, bed roughness and rheology.

2 Information within the correlation score surface

The backbone of velocity extraction from imaging satellites is image correlation (a.k.a. pattern matching and feature tracking). For a general overview of an image-matching pipeline, see Appendix A. The implementation of image correlation is done through the use of a subdomain or kernel in one image that is compared against a second image to find the most similar signal within this subdomain. Typically, a matching domain is a two-dimensional space (i,j), where each axis describes one translation. This leads to a two-dimensional correlation surface of similarity scores (Θ_i,j), where the highest score is taken as the candidate for the displacement.

Apart from the displacement information, other metrics can also be extracted from the correlation surface. For an extensive assessment of such metrics, see Xue et al. (2014). We interpret these not as metrics for dispersion but describe other qualitative aspects. For example, we interpret the absolute value of the highest peak as a proxy for confidence, while an indication for validity can be calculated from the ratio between the highest and second highest peak. Similarly, the signal-to-noise ratio is a proxy for uniqueness, but neither of these descriptors give any information about the matching precision. The just-mentioned reliability proxies are typically provided on a point-per-point basis within glacier velocity products, while an individual dispersion estimate is still lacking.

However, upon close inspection the width and form of the highest peak in a correlation surface changes and depends greatly on the image structure, for example, surfaces with a preferred orientation, such as crevasses (Fig. 1a). Here the maximum score is situated on a ridge of similar high scores, as there is a lack of distinguishable features along the direction of the crevasses. In the direction perpendicular to the dominant feature orientation, the correlation peak is sharp with steep flanks. In the direction of the feature orientation, though, the peak is weakly defined in one of the two directions and thus uncertain. Such correlation ridges occur abundantly on glaciers, as elongated features such as crevasses, moraines and streamlines populate many glacier surfaces. Paradoxically, it is these features that exhibit high contrast and are persistent over time and thus represent the dominant features for glacier displacement estimation.

Figure 1 Examples of cross-correlation results with anisotropy, due to the image content or the underlying surface flow. For both panels, the lower-left panel is the image template to be correlated with the upper-left search area. The resulting correlation surface is displayed in the lower right, and a zoom around the correlation peak is in the upper right. Both examples (a) and (b) show typical glacier surfaces that result in elongated correlation peaks (here called correlation ridges), sharply defined in one direction but weakly in the perpendicular direction. NCC: normalized cross-correlation.

A second process influencing the spread of correlation scores is when significant shear occurs within the template (Fig. 1b). Even though the variation (or contrast) in the template might be present in all directions, a ridge in the correlation surface can emerge similarly to the one from crevasses. In case of glacier surface shear, the simple translation assumed in the pattern matching is not valid (Debella-Gilo and Kääb, 2012). Misalignment in the outer parts of the template causes dissimilarity so that the correlation peak gets dampened and neighboring values increase at the same time, weakening the relative strength of the peak. If the size of the template is reduced in such situations, this spreading of correlation scores is reduced, however at the cost of a decreasing signal-to-noise ratio. A second remedy to shear or rotation is to apposition an affine model instead of one with translation only. Though non-linear and thus iterative in nature, such a higher-order model creates the opportunity to estimate shear and strain rates directly from the image matching (Debella-Gilo and Kääb, 2012).

Formulating the precision of a match can be done by looking directly at the variation in intensities within an image (Kanazawa and Kanatani, 2003). Local derivative filters can be used to describe the spatial variation within a Hessian matrix. However at which scale these filters should be set is in the case of naive image matching not always known. Nonetheless an image-based approach for precision estimation is beneficial when the scale is known, which is the case for feature descriptors such as scale-invariant feature transform (SIFT) or speeded-up robust features (SURF), and such formulations have been worked out (Zeisl et al., 2009). Another approach, which is also taken in this study, is to directly look at the correlation peak. Similarly to the image-based case, the curvature of the correlation peak can be described by the Hessian approach. This approach is implemented in Ampcorr an SAR-offset (synthetic-aperture radar) procedure within the ROI_PAC package (Repeat Orbit Interferometry Package; Rosen et al., 2004) and is described in a bit more detail in Casu et al. (2011). Here a similar approach is taken, but we directly model the correlation peak to a Gaussian function. From the background given in this section, our motivation is to capture information about the correlation surface, and in particular its peak, potentially allowing for a better judgment of the quality and precision of individual matches for displacement measurement.

3 Methodology

We perceive the close surroundings of the correlation function as a probability density function. This is a standard perception in the field of fluid mechanics (Bhattacharya et al., 2018), where pattern matching is known as particle image velocimetry (PIV). However, within this latter field of typically controlled laboratory environments, image matching is performed on small distinct features; hence shear effects are not present in the templates, while such effects are present for ice flow. Because the correlation surface is perceived as a probability density function, it is here fitted with a Gaussian function to be in line with generalized least-squares inversion techniques.

3.1 Co-variance from correlation spread

Here, we draw up a linear formulation to describe the variance of the correlation peak, which also considers its orientation. At a certain location in this search space (i,j), a two-dimensional Gaussian can be calculated through

$$\begin{array}{}\text{(1)}& f(i,j)=A\cdot e^{\left(-(a\cdot (i-i_{\mathrm{0}}{)}^{\mathrm{2}}+\mathrm{2}b\cdot (i-i_{\mathrm{0}})\cdot (j-j_{\mathrm{0}})+c\cdot (j-j_{\mathrm{0}}{)}^{\mathrm{2}})\right)},\end{array}$$

where i₀ and j₀ denote the center of the peak ($i_{\mathrm{0}}=i_{\max }+\mathrm{\Delta }i$) which might not coincide with the integer-valued location of the highest value in the correlation grid (i_max). The center of the top can be estimated by a peak-finding function and is here considered to be known. Equation (1) is in a simplified form, where A encompasses the magnitude and a, b and c are lumped constants. A detailed derivation thereof is given in Appendix B. Then, the rest of the unknowns can be estimated after some rearrangement:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \underset{\mathbf{y}}{\underbrace{\left[\begin{array}{c}\ln {\mathrm{\Theta }}_{i-\mathrm{1},j-\mathrm{1}}\\ \ln {\mathrm{\Theta }}_{i-\mathrm{1},j}\\ \mathrm{⋮}\\ \ln {\mathrm{\Theta }}_{i+\mathrm{1},j+\mathrm{1}}\end{array}\right]}}=\\ & \underset{\mathbf{A}}{\underbrace{\left[\begin{array}{cccc}\mathrm{1}& (i-\mathrm{1}-i_{\mathrm{0}}{)}^{\mathrm{2}}& (i-\mathrm{1}-i_{\mathrm{0}})\cdot (j-\mathrm{1}-j_{\mathrm{0}})& (j-\mathrm{1}-j_{\mathrm{0}}{)}^{\mathrm{2}}\\ \mathrm{1}& (i-\mathrm{1}-i_{\mathrm{0}}{)}^{\mathrm{2}}& (i-\mathrm{1}i{-}_{\mathrm{0}})\cdot (j-j_{\mathrm{0}})& (j-j_{\mathrm{0}}{)}^{\mathrm{2}}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \\ \mathrm{1}& (i+\mathrm{1}-i_{\mathrm{0}}{)}^{\mathrm{2}}& (i+\mathrm{1}-i_{\mathrm{0}})\cdot (j+\mathrm{1}-j_{\mathrm{0}})& (j+\mathrm{1}-j_{\mathrm{0}}{)}^{\mathrm{2}}\end{array}\right]}}\\ & \underset{\mathbf{x}}{\underbrace{\left[\begin{array}{c}\ln A\\ a\\ \mathrm{2}b\\ c\end{array}\right]}}.\end{array}\end{array}$$

This gives the possibility for directly estimating the unknowns (in x) through least-squares adjustment from the similarity scores (Θ). The direct neighborhood is used here (radius =1) when the peak is next to the border; otherwise a two-pixel radius is used. The lumped constants (a, b and c) can then be reformulated to extract the variances (σ²) and their dependency (ρ) from Eq. (1):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{2}\mathit{\rho }& ={\displaystyle \frac{b}{\sqrt{a\cdot c}}},{\mathit{\sigma }}_i^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{-\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot a}},\\ & {\mathit{\sigma }}_j^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{-\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot c}}.\end{array}\end{array}$$

This estimation procedure is an extension of Anthony and Granick (2009), who only resolved for σ_i and σ_j. However in the formulation of Eq. (1) the axes can have a dependency (ρ), and correlation ridges with different orientations can thus be estimated. Then the dispersion matrix (Q_yy) is composed of the estimates from Eq. (3) and the pixel spacing (d) as follows:

$$\begin{array}{}\text{(4)}& {\mathbf{Q}}_{yy}=\left[\begin{array}{cc}{d}_x& \mathrm{0}\\ \mathrm{0}& d_y\end{array}\right]\cdot \left[\begin{array}{cc}{\mathit{\sigma }}_i^{\mathrm{2}}& \mathit{\rho }{\mathit{\sigma }}_i{\mathit{\sigma }}_j\\ \mathit{\rho }{\mathit{\sigma }}_i{\mathit{\sigma }}_j& {\mathit{\sigma }}_j^{\mathrm{2}}\end{array}\right].\end{array}$$

The dispersion matrix (Q_yy) can directly be inserted into a co-variance matrix for error propagation or data assimilation. The off-diagonal elements of this matrix describe the dependencies between observations. Typically these are set to zero for displacement couples (e.g., Derkacheva et al., 2020), but they have the ability to describe the temporal and/or spatial relational dependencies within the dataset (a.k.a. spatial coherency; Riel et al., 2014).

3.2 From (co-)variance to standard-error ellipse

For the dependencies between two-dimensional displacements, as presented here, interpretation of the elements within the dispersion matrix might not be intuitive. For example, an equal variance can still produce an orientation dependency, as can be seen for example in Fig. B1. Hence, here we give the transformation from the standard-error axis (${\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}},{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}$) to a description of standard-error ellipse in the form of a minor and a major axis (λ₁,λ₂, respectively) and its orientation (θ). The two axes can be extracted through

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{2}}}+\sqrt{{\displaystyle \frac{({\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}-{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}{)}^{\mathrm{2}}}{\mathrm{4}}}+\mathit{\rho }{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}},\\ & {\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{2}}}-\sqrt{{\displaystyle \frac{({\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}-{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}{)}^{\mathrm{2}}}{\mathrm{4}}}+\mathit{\rho }{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}.\end{array}\end{array}$$

Similarly, the orientation of the ellipse (θ) can be calculated by

$$\begin{array}{}\text{(6)}& \tan \left(\mathrm{2}\mathit{\theta }\right)={\displaystyle \frac{\mathrm{2}\mathit{\rho }{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}-{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}.\end{array}$$

3.3 Derivatives of flow from incomplete data

Surface strain rates are used in this study to assess the relation between the correlation ridge and ice deformation. Such strain rates can be extracted from a velocity field; however remote sensing results contain holes and patches without estimates, since similarity could not be established. Hence a robust estimation framework is given in Appendix C that is somewhat resistant to such sporadic outliers. This procedure is used here to have a more complete strain rate field for analysis.

3.4 Crevasse characteristics from a Radon transform

To assess the impact of directionality in the input images on our approach to compute and use the dispersion of individual correlations, we need to quantify the directional characteristics of glacier images. In particular crevasse fields have strong directional properties, which can be composed of cracks with several predominant orientations. In order to extract the local crevasse characteristic for each matching template, a Radon transform is used, as described in earlier work (Gong et al., 2018). This methodology provides an argument of the strongest crevasse direction and a strength of this signal. With both the shear flow and crevasse orientation quantified, it is then possible to assess the sensitivity of image matching to these two properties.

Figure 2 Color-coded speed and streamlines of Sermeq Kujalleq (Jakobshavn Isbræ) between 20 and 30 July 2020 based on Sentinel-2A imagery. The upper inset (b) shows the along-flow variance, while the lower inset (c) shows the across-flow variance.

4 Results

Here we present results from two sites, namely Sermeq Kujalleq (Jakobshavn Isbræ), Greenland, and Malaspina Glacier system, Alaska.

4.1 Sermeq Kujalleq, Greenland

We demonstrate and assess our method to estimate the uncertainty in displacement matching using a small subset of two orthorectified Sentinel-2A scenes over Sermeq Kujalleq, a large and fast outlet glacier of the Greenland ice sheet. High-pass filtered imagery (following Fahnestock et al., 2016) of the 10 m near-infrared band number 4 is used. The image pair has acquisitions that are 10 d apart, acquired from the same orbit. We apply a template window of 200 m in dimension, and velocities are estimated every 100 m, with a search window of 800 m. Co-registration is not applied to the image pair beforehand, as related offsets are not in the scope of this study nor does the absence of co-registration influence the outcomes of the presented work. A similar template size of 200 m was used for the crevasse detection using the Radon transform.

The velocity magnitude between the two images (20 and 30 July 2020), derived streamlines, and the resulting along- and across-flow variance estimates are shown in Fig. 2. The streamlines indicate a strongly convergent flow of this outlet glacier. At some places the signal-to-noise ratio of the image matching was too low (SNR<4), and such displacements have been excluded. This happened in particular at the eastern part, where a cloud is present in one of the images, and at other locations which seem to correspond to supraglacial lakes. Along-flow variations are large at the northern side of the main outlet and in the bend before the outlet terminates in the fjord. Across-flow variation occurs in the more slowly moving regions, where sheared crevassing occur, such as the terminus of the outlet northwest and the southwestern part of the study region.

Figure 3 Descriptions of directionality for the case study of Sermeq Kujalleq. (a) The correlation peak orientation for individual image-matching results. (b) The direction of the imagery, extracted from the Radon transform, that directly operates upon the imagery.

The dominant crevasse orientation (Fig. 3b) is transverse to the flow direction, aligning with crevasses originating from extensive extensional strain. This has a stark similarity with the orientation of the correlation surface (Fig. 3a). Some regions have more complex orientations, most likely due to variations in surface slope and bedrock variation.

Figure 4b shows the shear strain rate. A major feature is a large shear zone along the southern flank of the main flow channel. Finer details like alternating patches are also present in the main outlet, which could stem from the propagation of subglacial features to the surface. In Fig. 4a, a measure for the elongation of the correlation ridge is plotted. Here, elongation is given as the normalized inequality of the two dispersion components ($[min({\mathit{\lambda }}_x,{\mathit{\lambda }}_y)-max({\mathit{\lambda }}_x,{\mathit{\lambda }}_y\left)\right]/[{\mathit{\lambda }}_x+{\mathit{\lambda }}_y]$). Hence 0 corresponds to a perfectly circular distribution, while 1 would be a straight ridge.

Figure 4 The strength of asymmetry of the correlation peak (a) and estimated surface shear (b).

4.2 Malaspina Glacier system, Alaska

Results from the surroundings of Malaspina Glacier and Agassiz Glacier in the St. Elias Mountains are presented here as well. The region exhibits more supraglacial features than Sermeq Kujalleq, which is an outlet of the Greenland ice sheet with predominantly clean ice. For example, a large collection of moraines, ogives, foliations, meltwater channels and more diverse orientations of flow is present on both Malaspina Glacier and Agassiz Glacier, as can be seen in Fig. D1 in the Appendix. Malaspina Glacier is an outlet of Seward Glacier with a total area of 5000 km² (Molnia, 2008); its entire piedmont lobe lies within the ablation area. Agassiz Glacier is the other large tributary of the Malaspina Glacier system and creates a distinct western lobe. The ice transport from Seward Glacier has multi-annual fluctuations (see https://www.youtube.com/watch?v=YslhQZwvvu0, last access: 1 June 2022, and Altena et al., 2019), creating looped or curved moraine bands.

Here we use two subsets of Sentinel-2 scenes from 21 August and the 15 September 2019, a 25 d difference, and from the same orbit. Processing parameters are similar to the Sermeq Kujalleq study: a high-pass-filtered band 4 image was matched, with a template window of 200 m wide, and velocities are estimated every 100 m, with a search window of 800 m. No co-registration over stable ground was done, so the velocities should be seen as displacement (being real surface displacements or artificially created due to sensor/processing biases).

Figure 5 Image displacement (a) and elongation of the correlation ridge (b). The red parallelogram illustrates large values of correlation spread, which also aligns with shear flow as in Fig. 6a, while the large spread highlighted by the red circle is probably due to crevassing, since flow is fairly homogenous in this region.

The estimated displacements over the study region (Fig. 5a) have a smooth surface. In the mountains region, small speckles are present, as well as a small patch on Agassiz Glacier. In this zone the transient snow line was located, so correspondence is more difficult to establish. Similarly, haze or thin cloud cover might be present at the end of the snout of Agassiz Glacier, causing further correspondence failures.

Figure 6 Estimated surface shear, derived from the estimated velocity (Fig. 5a). The red parallelogram illustrates a region with clear shear flow, but also other regions of Malaspina Glacier have large shear rates and correspond to large elongations. Red arrows within the large circle in the center highlight faint shear patterns that stem from sensor misalignment (Fig. D2b). The strength of image structure in the form of the Radon transform (b) has some distinct features of interest; the red rectangle highlights strong crevasses, which also result in elongated correlation spread (Fig. 5b). The red circle at Agassiz Glacier highlights another region with strong crevasses, which also is present in the estimates of the flow divergence (Fig. 10a) and the signal-to-noise estimate of the image matching (Fig. 9b).

The pattern of elongation (Fig. 5b) and shear (Fig. 6a) are similar at the borders of Agassiz Glacier, as indicated by the red parallelogram. This is not the case for many other parts, while a region which shows no extensive local shear or extension can be seen at the start of the lobe (encircled in red). However, this region does have heavy crevassing as well (Fig. 6b). This is not only happening in this region, but in general the Radon strength correlates well with the elongation of the ridge. Hence, excessive shear and extension might create crevasses, and these seem to be the most dominant mechanism for asymmetrical correlation spread. Other signals are also present in the shear estimate (Fig. 6a), but these will be highlighted later in the discussion, as they are not related to correlation spread.

4.2.1 Orientation of crevasses and dispersion peak

The dominance of the feature orientation (Fig. 7a) to the direction of the correlation ridge (Fig. 7b) is present here, as is also observed in the Sermeq Kujalleq case. The structure of these two independent proxies are very similar. While Sermeq Kujalleq is dominated by clean ice, it seems also other directional features like foliations and moraine bands influence the correlation surface.

Figure 7 Orientation descriptors over the Malaspina case study, estimated through Radon transform (Fig. 9a) and correlation spread (Fig. 9b).

5 Discussion

5.1 Interpretation of the dispersion signal

In general the main orientation of the crevasses at Sermeq Kujalleq (Fig. 3b) seems to correspond to the orientation of the Gaussian peak (Fig. 3a). When these two parameters are plotted against each other their relation becomes even clearer (Fig. 8a). The bulk of crevasse orientations are oriented towards a north–south axis, corresponding to be perpendicular to the main flow direction. A straight correlation between both parameters is present in Fig. 8a but does not cover the whole domain equally due to the limited distribution of flow directions. A relation with crevasse presence is profound (Fig. 8b); when the Gaussian peak is close to symmetrical (i.e., inequality near zero), there is no clear relation, but this increases when crevasses become more apparent in the imagery through the Radon transform. The pattern of elongation of Sermeq Kujalleq (Fig. 4a) is less pronounced and does not have a clear linear relation to shear flow (not shown). A reason why no clear relation between shear flow and elongation of the correlation peak is found in our example can be due to the strong presence of crevasses that then dominate the signal and image correlation in this dataset. The dominance of crevasses in the study region could suppress the existence a clear relation with non-uniform ice flow. Nonetheless, crevasses seem to be the dominant driver for asymmetry in the correlation peak.

Figure 8 Probability density plots of results for Sermeq Kujalleq of correlation peak versus crevasse orientation (Figs. 2 and 3a) and asymmetry of the correlation peak and versus crevasse strength (Figs. 2 and 4a).

Figure 9 Correlation descriptors over the Malaspina case study showing the absolute correlation value for each match (a) and the signal-to-noise values (b). The red circles indicate highly distinct regions that align with regions of large shear (Fig. 6a) and many crevasses (Fig. 6b); the red square also indicates a region with high values but has homogenous flow with many surface cracks. The red diamond highlights a pattern that is similar to the integer distance as shown in Fig. 10b, indicating a dependency.

5.2 Description of dispersion

In earlier work the handing of dispersion has been estimated through sampling statistics (standard deviation and mean absolute difference), where displacement estimates are compared against in situ measurements or stable terrain. The use of stable terrain for dispersion estimation has drawbacks, apart from assuming constant variance of the whole scene as mentioned earlier. Specifically, image matching in the frequency domain is hampered by peak locking that favors integer displacements (Foroosh et al., 2002); thus in the configuration of stable terrain (in this case zero displacement), sample statistics will give an opportunistic estimate of precision. A dispersion formulation based on image intensities has been proposed (Förstner, 1987). Then template matching itself is formulated within a least-squares framework where the noise level of individual pixels propagates into a precision estimate of a match, but such estimates can be highly influenced by sample statistics where a large amount of pixels in a template cause the system of equations to produce a very good measurement precision; furthermore outliers in such a formulation are neglected (Maas et al., 2010).

Thus the method presented here can be a direction to formulate measurement precision, without biases introduced by sample statistics and peak locking. Another advantage of our method is the possibility for using statistical testing (Teunissen, 2000) and integration into data assimilation models or time series construction through a richer description of the co-variances (Riel et al., 2021).

We postulate that the correlation coefficient is a proxy for the confidence of a match and are therefore less suited to function as a descriptor of precision. The maximum correlation coefficient and the signal-to-noise proxy are dissimilar proxies. For example, the narrow and crevassed outlet of Malaspina Glacier has low correlation scores (Fig. 9a) but a high signal-to-noise ratio (Fig. 9b). Upon closer inspection, a striking feature of multiple peaks grouped together might be observed in the correlation score (see the region inside the red diamond). This pattern aligns with the sub-pixel displacement away from an integer, as the correlation score is estimated at individual steps. This off-integer bias in the correlation score can be replicated by the displacement estimate and is done so in Fig. 10b. Hence, using a correlation score as precision proxy (e.g., Ding et al., 2021), while it is contaminated by off-integer biases, is not recommended.

A second commonly used proxy for precision is the signal-to-noise ratio. Here we postulate that this proxy might describe the uniqueness of a match. Very high signal-to-noise values (Fig. 9b) seem to coincide with strong crevassing (Fig. 6b), as is also indicated by the red encircling. A second class of high values (see red square) is present in clean ice zones of the lobe of Malaspina Glacier, where distinct foliations occur, giving a unique fingerprint for the matching.

Figure 10 Estimated surface divergence, derived from the estimated velocity (Fig. 5a). The red circles indicate regions where significant crevassing is present (Fig. 6b). The modulus from a combination of sub-pixel displacements (Fig. D2a, b) is shown in (b). The red diamond indicates a pattern that is similar to the absolute correlation value as shown in Fig. 9a.

In this study we propose to use a Gaussian formulation to describe the matching precision. If the maximum correlation or signal-to-noise ratio would be a good proxy for precision, then one can expect a correlation or some form of agreement between the major axis (Fig. 11) and these other proxies. However, for the data over Malaspina Glacier this does not seem to be the case, as these proxies only seem to be correlated in the extreme ends. Thus, the proposed dispersion parameters do provide a new type of data description, which we think has a straightforward connection to measurement precision.

Figure 11 Probability scatterplots between different matching descriptors for the Malaspina Glacier system.

5.3 Implementation issues

The implementation done here for our correlation-dispersion-based method is a simple least-squares adjustment, and no robust re-weighting is applied. This can result in negative variances or rank deficiency, corresponding to the white data voids in Fig. 3. Causes for such anomalies can come from the logarithmic function within Eq. (2), capable of transforming white to asymmetric noise (Anthony and Granick, 2009).

In this study, the correlation computation is done in the spatial domain. When transformed back to the spatial domain, frequency domain methods produce sharp peaks in the correlation surface in the form of a two-dimensional Dirichlet function, as they prescribe consistent rigid displacement at integer resolution (Foroosh et al., 2002). Furthermore, when sufficient shear occurs or repeating image features are present, this might result in multiple distinct but sharp peaks in the correlation surface (Scarano, 2001). Hence, interpretation of our dispersion estimation is most suited for spatial-domain methods.

Finally to demonstrate its application domain, we introduce a generalized least-squares framework to use our dispersion estimation (see Appendix C) and resolve issues caused by missing data from neighboring displacement estimates when estimating strain rates. This is a step towards a more integrated approach and moves away from parameter-based interpolation (e.g., Lüttig et al., 2017).

6 Conclusions

Quantifying the measurement precision of individual displacement estimates from matching repeat spaceborne images has received little attention in recent years despite the increasing efforts to produce large displacement datasets from an increasing number of suitable data. Here, we introduce a simple procedure to estimate the correlation dispersion of such displacement measurements (either optical or SAR), through characterizing the shape of the correlation surface. We demonstrate this technique for Sermeq Kujalleq, a fast-flowing and heavily crevassed outlet of the Greenland ice sheet and the Malaspina Glacier system. Dispersion results are compared to shear strain rates and crevasse orientation. These results indicate that crevasses are the dominant driver for asymmetry in the correlation surface. We suggest this simple procedure to estimate uncertainty in individual image matches can be useful in processing pipelines for large-volume image displacement measurements, so error-propagation can be applied on a large scale and will improve inversion of other geophysical properties. In all, we hope this demonstrates the rich information present in satellite imagery and its processing chain and might make it easier to extract a more detailed physical signal from such noisy remote sensing products.

Appendix A: Schematic of an image-matching pipeline

In order to clarify where in a displacement processing scheme the dispersion estimation can be implemented, a schematic of an image-matching pipeline is drawn in Fig. A1. It consists of the following five steps:

Figure A1 Schematic of the main procedure to generate a displacement field from a pair of remote sensing images.

• i.

  Given the extent of the imagery, a mask is generated indicating what is ocean, land and glacier.

• ii.

  A regular grid is generated, where for each location the land cover is recorded.

• iii.

  For each post of the grid, a subset of the satellite imagery is used. A kernel is moved over a base image, and at every location a similarity score is estimated. This generates a correlation surface. The highest value is taken as the correct displacement. The neighboring correlation values of this peak can be used for sub-pixel localization, but the same values can also be used for the dispersion calculation following the method presented in this study.

• iv.

  The displacements over stable ground are used to correct offsets due to misalignment of the satellite platform.

• v.

  The co-registration parameters are subtracted from the displacement vectors, resulting in a grid of velocities and its precision.

Appendix B: Complete derivation

A two-dimensional normal distribution, with a dependency (ρ) in its variables can be written as (Teunissen et al., 2009)

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}I& (x,y)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_x{\mathit{\sigma }}_y\sqrt{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}}\cdot \exp \left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\right.\\ & \left.\left({\displaystyle \frac{(x-x_{\mathrm{0}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}_x^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{2}\mathit{\rho }(x-x_{\mathrm{0}})(y-y_{\mathrm{0}})}{{\mathit{\sigma }}_x{\mathit{\sigma }}_y}}+{\displaystyle \frac{(y-y_{\mathrm{0}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}}}\right)\right],\end{array}\end{array}$$

where x and y denote coordinates on two orthogonal axis, σ² is the variance, and x₀ and y₀ are their mean. This formulation can be written out fully with parameters (A, a, b and c) substituted for the ease of readability (Eq. 1). This results in a linear system of equations with four unknowns, so these need to be estimated through several neighboring correlation values, as written down in Eq. (2). The following operations show the transformation from one formulation to the other.

$$\begin{array}{}\text{(B2)}& {\displaystyle }\begin{array}{rl}I(x,y)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_x{\mathit{\sigma }}_y\sqrt{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}}\exp \left[\left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot {\mathit{\sigma }}_x^{\mathrm{2}}}}\right.\right.\\ & \cdot (x-x_{\mathrm{0}}{)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{2}\mathit{\rho }}{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot {\mathit{\sigma }}_x{\mathit{\sigma }}_y}}\cdot (x-x_{\mathrm{0}}\left)\right(y-y_{\mathrm{0}})\\ & \left.\left.-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot {\mathit{\sigma }}_y^{\mathrm{2}}}}\cdot (y-y_{\mathrm{0}}{)}^{\mathrm{2}}\right)\right]\end{array}\text{(B3)}& {\displaystyle }\begin{array}{rl}I(x,y)& =A\cdot \exp \left[a\cdot (x-x_{\mathrm{0}}{)}^{\mathrm{2}}+b\cdot (x-x_{\mathrm{0}})\right.\\ & \left.\cdot (y-y_{\mathrm{0}})+c\cdot (y-y_{\mathrm{0}}{)}^{\mathrm{2}}\right]\end{array}\text{(B4)}& {\displaystyle }\begin{array}{rl}\ln \left[I(x,y)\right]& =\ln \left[A\right]\cdot \mathrm{1}+a\cdot (x-x_{\mathrm{0}}{)}^{\mathrm{2}}+b\cdot (x-x_{\mathrm{0}})\\ & \cdot (y-y_{\mathrm{0}})+c\cdot (y-y_{\mathrm{0}}{)}^{\mathrm{2}}\end{array}\end{array}$$

The substituted parameters (A, a, b and c) can be written out fully as

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}A& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_x{\mathit{\sigma }}_y\sqrt{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}},a=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot {\mathit{\sigma }}_x^{\mathrm{2}}}},\\ & b={\displaystyle \frac{\mathrm{2}\mathit{\rho }}{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot {\mathit{\sigma }}_x{\mathit{\sigma }}_y}},c=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot {\mathit{\sigma }}_y^{\mathrm{2}}}}.\end{array}\end{array}$$

Transferring these lumped parameters towards the Gaussian parameters (Eq. B1) is done though first formulating them in relation to the dependency (ρ):

$$\begin{array}{}\text{(B6)}& {\displaystyle }\mathrm{2}\mathit{\rho }={\displaystyle \frac{b}{\sqrt{a\cdot c}}},\text{(B7)}& {\displaystyle }\begin{array}{rl}a\cdot c& =-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot {\mathit{\sigma }}_x^{\mathrm{2}}}}\cdot -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot {\mathit{\sigma }}_y^{\mathrm{2}}}}\\ & ={\displaystyle \frac{\mathrm{1}}{{\mathrm{2}}^{\mathrm{2}}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}{)}^{\mathrm{2}}\cdot {\mathit{\sigma }}_x^{\mathrm{2}}\cdot {\mathit{\sigma }}_y^{\mathrm{2}}}},\end{array}\text{(B8)}& {\displaystyle }\begin{array}{rl}\sqrt{a\cdot c}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot {\mathit{\sigma }}_x\cdot {\mathit{\sigma }}_y}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot {\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot {\mathit{\sigma }}_x\cdot {\mathit{\sigma }}_y}},\\ & b={\displaystyle \frac{\mathrm{2}\mathit{\rho }}{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot {\mathit{\sigma }}_x\cdot {\mathit{\sigma }}_y}}.\end{array}\end{array}$$

Knowing ρ makes it possible to solve the other equations and extract the variances (${\mathit{\sigma }}_x^{\mathrm{2}}$, ${\mathit{\sigma }}_y^{\mathrm{2}}$) from the other lumped parameters (Eq. 3):

$$\begin{array}{}\text{(B9)}& \begin{array}{rl}{\mathit{\sigma }}_x^{\mathrm{2}}& ={\displaystyle \frac{\mathrm{1}}{-\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot a}}={\displaystyle \frac{-\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}{-\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\cdot {\mathit{\sigma }}_x^{\mathrm{2}},\\ {\mathit{\sigma }}_y^{\mathrm{2}}& ={\displaystyle \frac{\mathrm{1}}{-\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\cdot c}}={\displaystyle \frac{-\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}{-\mathrm{2}\cdot (\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\cdot {\mathit{\sigma }}_y^{\mathrm{2}}.\end{array}\end{array}$$

With the resulting parameters (σ₁, σ₂ and ρ), an oriented ellipse can be described, as shown in Fig. B1.

Figure B1 Example of ellipses with different dispersion parameters. Illustration adopted from Polman and Salzmann (1996).

Appendix C: Kernel computation in a generalized least-squares framework

Flow descriptors like strain rates can also give an insight into the geometric bedrock configuration or properties related to subglacial sliding. Strain rates can be formulated in relation to the local flow direction, giving longitudinal, transverse or shear flow, respectively. These properties are computed from velocity estimates over a close neighborhood of surrounding pixels. As strain rates are derivatives of velocities, they are particularly sensitive to the propagation of noise and errors in the input velocities. Applying thresholds and filters to the strain rates based on variations or low quality of the input velocities can lead to voids in the resulting strain rate field. Here, a methodology is introduced that is somewhat resistant to such cases caused by velocity errors or missing data.

The methodology presented here is based upon the redundancy of a kernel, since it is typically formulated as a smoothed differentiation. The steps are schematically illustrated in Fig. C1. A convolution (⊗) is written out as a matrix form of a displacement grid (P) and a kernel (G). In this matrix form one can see it clearly as a weighted linear combination from neighboring velocity measurements. For the sake of clarity, the examples shown in this schematic are an implementation of two different kernels (a Sobel and Prewitt), for the two different spatial axis (x, y). Each column in the design matrix (A) is independent and is composed of positive and negative entities. The summation of all elements within the kernel need to cancel each other out, as is indicated by the colored elements. However, when gaps occur in the neighborhood, this energy balance is disrupted. Consequently, this lost weight should be added to others within its group or, reversely, taken away from entries with the group with an opposing sign. When the convolution is written out directly in matrix form, this allocation of energy is done by column-wise operators. If the neighborhood is out of balance, the kernel is not estimated.

Figure C1 Schematic of a computation of a convolution, in this case the first derivative in the vertical and horizontal direction.

In the example shown in Fig. C1 the horizontal and vertical components (x, y) are independent. However the dependency can also be included, since formulating a convolution as a least-squares estimation makes it possible to propagate the co-variances. Hence, the co-variances of the image matching as given in Eq. (4) can be used to estimate the precision and dependencies of derived parameters, through

$$\begin{array}{}\text{(C1)}& {\mathbf{Q}}_{\widehat{x}\widehat{x}}={\mathbf{A}}^{\top }{\mathbf{Q}}_{yy}^{-\mathrm{1}}\mathbf{A},\end{array}$$

hence, estimating derivatives with correct weighting, making generalized least squares possible as

$$\begin{array}{}\text{(C2)}& \widehat{\underset{\mathrm{¯}}{\mathbf{x}}}={\mathbf{Q}}_{\widehat{x}\widehat{x}}^{-\mathrm{1}}{\mathbf{A}}^{\top }{\mathbf{Q}}_{yy}^{-\mathrm{1}}\underset{\mathrm{¯}}{\mathbf{y}}.\end{array}$$

Nevertheless, improvement is only made on a local level in a direct neighborhood covered by the kernel, so when large parts are affected with regions of missing values or when the outlier detection is false, spurious fluctuations can still propagate into the final product.

Appendix D: Additional information about the Malaspina Glacier case

In this appendix, additional illustrations are shown for the Malaspina Glacier to ease interpretation of the results and to highlight the information present in dispersion peak.

Here we also show some sub-pixel displacement plots, as the integer (Fig. 10b) is based upon the combination of two axes, as the remainder of the modulus of displacement is shown in Fig. D2. Sensor specific artifacts are present in these figures as indicated by the gray boxes or the red encircled region; see Kääb et al. (2016) and Stumpf et al. (2018) for more details. They are mentioned here explicitly, since these patterns are striking and might dilute the interpretation of the results in Fig. 6.

Figure D1 Sentinel-2 scene over Malaspina Glacier (center) and Agassiz Glacier (left) with annotations in red to enhance interpretation.

Figure D2 Rainbow-color-coded remainder of the modulus of displacement, for the horizontal and vertical direction (a, b). The gray arrows and boxes show displacement artifacts due to internal sensor alignments. The red circles and arrows indicate oscillations that stem from similar internal registration errors.