the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Coregistration and residual correction of digital elevation models: a comparative study
Tao Li
Yuanlin Hu
Bin Liu
Liming Jiang
Hansheng Wang
Xiang Shen
Digital elevation models (DEMs) are currently one of the most widely used data sources in glacier thickness change research, due to the high spatial resolution and continuous coverage. However, raw DEM data are often misaligned with each other, due to georeferencing errors, and a coregistration procedure is required before DEM differencing. In this paper, we present a comparative analysis of the two classical coregistration methods proposed by Nuth and Kääb (2011) and Rosenholm and Torlegard (1988). The former is currently the most commonly used method in glacial studies, while the latter is a seminal work in the photogrammetric field that has not been extensively investigated by the cryosphere community. Furthermore, we also present a new residual correction method using a generalized additive model (GAM) to eliminate the remaining systematic errors in DEM coregistration results. The performance of the two DEM coregistration methods and three residual correction algorithms (the GAMbased method together with two parametricmodelbased methods) was evaluated using multiple DEM pairs from the Greenland Ice Sheet and mountain glaciers, including Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEMs, ZiYuan3 (ZY3) DEMs, the Shuttle Radar Topography Mission (SRTM) DEM, and the Copernicus DEM. The experimental results confirm our theoretical analysis of the two coregistration methods. The method of Rosenholm and Torlegard has a greater ability to remove DEM misalignments (an average of 4.6 % and 13.7 % for the test datasets from Greenland Ice Sheet and High Mountain Asia, respectively) because it models the translation, scale, and rotationinduced biases, while the method of Nuth and Kääb considers translation only. The proposed GAMbased method performs statistically better than the two residual correction methods based on parametric regression models (highorder polynomials and the sum of the sinusoidal functions). A visual inspection reveals that the GAMbased method, as a nonparametric regression technique, can capture complex systematic errors in the DEM coregistration residuals.
 Article
(16744 KB)  Fulltext XML

Supplement
(6003 KB)  BibTeX
 EndNote
Differencing between multitemporal digital elevation models (DEMs) is a widely used approach for mapping glacier elevation changes at local and regional scales (Bolch et al., 2011; Gardelle et al., 2013; Pieczonka et al., 2013; Liu et al., 2019; Hugonnet et al., 2021). However, limited by the imaging and georeferencing techniques, systematic errors often exist in the raw DEMs (Rodriguez et al., 2006), which can lead to wrong estimation of glacier mass change and false detection of glacier surges (Nuth and Kääb, 2011). Numerous studies have confirmed that a coregistration process is required to remove these biases before DEM differencing is conducted (Van Niel et al., 2008; Nuth and Kääb, 2011; Paul et al., 2015).
DEM coregistration has been extensively studied, and the existing methods can be broadly classified into two main categories. The first category requires an explicit datamatching process (i.e., correspondence search). Typical methods in this category include featurepointbased methods, e.g., the scaleinvariant feature transform (SIFT) descriptor (Aguilar et al., 2012; Sedaghat and Naeini, 2018) and the method based on centroids of subwatersheds (Li et al., 2017); featurelinebased methods, e.g., methods based on stream networks or watershed boundaries (Karkee et al., 2008); the multifeaturebased surfacematching method (Wu et al., 2013); the iterative closest point (ICP) algorithm and its variants (Besl and Mckay, 1992; Rusinkiewicz and Levoy, 2001; Di et al., 2012); and the least squares 3D surface matching (LS3D) algorithm (Gruen and Akca, 2005; Akca, 2010). All of the above methods originate from image or point cloud processing studies, and they can be used for the coarse coregistration of DEMs without georeferenced information. However, the main disadvantage of these methods is that the correspondence finding procedure is very timeconsuming when processing large DEMs. Moreover, the accuracy of the imagebased methods (e.g., SIFT) is strongly dependent on extracting a large number of highquality features, which is not an easy task for DEMs lacking sufficient textures.
The second category of DEM coregistration methods does not require an explicit matching process. The optimization objective of these methods is usually to minimize the sum of the vertical distances between two DEMs, where each pixel in the secondary DEM implicitly corresponds to the same planimetric position in the reference DEM. These methods are not suitable for images lacking georeferenced information, but they are strongly recommended for highaccuracy applications where the DEMs have been georeferenced or coarsely coregistered. The typical algorithms in this category include grid search methods (Hofton et al., 2006; Berthier et al., 2007; Van Niel et al., 2008; Cucchiaro et al., 2020) and terraininformationbased methods (Gorokhovich and Voustianiouk, 2006; Peduzzi et al., 2010; Nuth and Kääb, 2011). The grid search methods search for the best alignment result by stepwise shifting the secondary DEM a little bit, alternatively along x and y directions in a predefined window (e.g., 5×5 pixels). However, these methods have been rarely used in the recent literature because their bruteforce search process comes with a huge computational cost. The terraininformationbased methods are derived from the analytical relationship between the elevation differences of the DEMs and terrainrelated information. The method proposed by Nuth and Kääb (2011) employs the terrain slope and aspect as explanatory variables in the regression model and is currently the most commonly used DEM coregistration algorithm in glacial studies (Vacaflor et al., 2022). In a much earlier study, Rosenholm and Torlegard (1988) developed an absolute orientation algorithm for stereo models based on terrain gradients. This method has been widely used for DEM coregistration in the photogrammetry field but, unfortunately, has rarely been considered in the cryosphere community.
The goal of this paper is twofold. The first goal is to reveal the connections and differences between the slope/aspectbased method of Nuth and Kääb (2011) (hereinafter referred to as NK) and the terraingradientbased DEM coregistration algorithm of Rosenholm and Torlegard (1988) (hereinafter referred to as RT), and the second goal is to present a nonparametric approach to remove the complex systematic errors in DEM coregistration results.
This section focuses on the analytical (i.e., the terraininformationbased) DEM coregistration methods only. We will demonstrate that the NK and RT methods are theoretically compatible. As the original algorithms in the works of NK and RT were presented in distinct forms, we will present detailed derivations of the equations used in their algorithms and variants.
2.1 The method of Nuth and Kääb
2.1.1 Standard version
The equations of the NK method are derived from the geometric relationship (see Fig. 1) of the elevation differences induced by the DEM shift with respect to the terrain slope (θ) and aspect (ψ) values. Firstly, we consider the special case where b=ψ (where b is the aspect of the shift vector); i.e., the translation is exactly along the terrain aspect direction. As shown in Fig. 1b, the induced elevation difference is given by
where a and c are the horizontal and vertical distances of the shift vector, respectively.
In a more general scenario, b≠ψ. As shown in Fig. 1a, the horizontal shift vector OE' is decomposed into OE and EE'. Since EE' is perpendicular to the vertical plane OEF defined by the gradient vector and the terrain aspect direction, this does not cause any elevation change. The vertical difference induced by OE' is therefore equal to that of OE, and we have
The above equation is Eq. (2) in Nuth and Kääb (2011). In this paper, we refer to this as the NK standard version. The cylindrical coordinates $\left(a,b,c\right)$ of the shift vector can be estimated from a nonlinear regression of Eq. (2), and the corresponding Cartesian coordinates are then given by
As shown in Fig. 2, an iterative process is generally required for accurate coregistration of two DEMs (Nuth and Kääb, 2011), where the coordinates of the secondary DEM are updated in every iteration by the following equation:
where the subscript i represents the ith iteration. The iterative process terminates when the change in the dispersion characteristics (median absolute deviation from zero) of the elevation differences between iterations is less than a predefined threshold.
2.1.2 Simplified version
Nuth and Kääb (2011) did not use Eq. (2) in their experiments but instead adopted a simplified regression equation by dropping one explanatory variable (θ). Firstly, both sides of Eq. (2) are divided by tan (θ):
θ in the right side of the equation is then approximately replaced by the mean terrain slope of the DEM:
where
Accordingly, the shift vector in the Cartesian coordinate system is given by
2.1.3 Linear version
The standard version of the NK method can be converted to a linear regression equation by combining Eq. (2) with Eq. (3):
This equation uses ΔX,ΔY, and ΔZ as the regression coefficients directly, and, accordingly, the conversion of the shift vector from cylindrical coordinates to Cartesian coordinates is no longer required.
2.2 The method of Rosenholm and Torlegard
In the RT method, the misalignment between two DEMs is described by a 3D similarity transformation (Kim and Jeong, 2011), and the coordinate update equation for the secondary DEM is
where γ is the scale factor; ω,φ, and κ are the rotation angles (in radians) about the X,Y, and Z axes, respectively; and the subscript C refers to the coordinates being zerocentered. Note that the values of $\mathit{\gamma},\mathit{\omega},\mathit{\phi}$, and κ are relatively small in the DEM coregistration process. The coordinate changes in each iteration can be approximated as
By comparing Eq. (11) with Eq. (4), it can be seen that the coordinate change of every DEM pixel is a constant vector ($\mathrm{\Delta}X,\mathrm{\Delta}Y,\mathrm{\Delta}Z$) in the NK method, while it varies with the position (${X}_{\mathrm{C}},{Y}_{\mathrm{C}},{Z}_{\mathrm{C}}$) in the RT method. Accordingly, the following equation is derived by substituting Eq. (11) into Eq. (9):
By rearranging the above equation, the DEM elevation differences caused by translation, scaling, and rotation can be obtained as
where
It can be found in the geoscience literature that the slope and aspect angles relate to the terrain gradients, and the following equation (Peckham and Jordan, 2007) exists when the terrain aspect is measured clockwise from north:
where f_{X} and f_{Y} are the gradients of the terrain in the X and Y directions, respectively. From the above equation, we obtain
Finally, by substituting Eq. (16) into Eq. (14), Eq. (13) is as follows:
The above equation is Eq. (6) in Rosenholm and Torlegard (1988). Comparing Eqs. (17) and (13) with Eq. (9), it can be seen that the RT and NK methods are theoretically compatible with each other, though the algorithms in their original works were presented in distinct forms.
A residual correction procedure is highly recommended after DEM coregistration (Berthier et al., 2007; Girod et al., 2017) because some systematic errors related to the terrain height and satellite acquisition geometry (alongtrack and crosstrack) often remain. As elevationdependent biases were not observed in our experiments, the following section introduces the residual correction algorithms for the alongtrack and crosstrack directions only.
3.1 Parametric regression
Highorder polynomial (sixth to eighth order) regression is the most commonly used way to fit DEM coregistration residuals (Nuth and Kääb, 2011; Gardelle et al., 2013; Berthier et al., 2016; Brun et al., 2017) and is usually performed in a stepwise manner:
with
where X_{t} and Y_{t} are the crosstrack and alongtrack coordinates, respectively; θ_{t} is the angle between the alongtrack direction and the north; m is the degree of the polynomial; and P_{i} and P_{j} are the coefficients to be estimated.
Previous studies have reported that the residual signals in the alongtrack direction often appear at one to three frequencies and are most likely induced by satellite attitude jitter, which is mainly caused by highfrequency mechanical vibration (Nuth and Kääb, 2011). Girod et al. (2017) pointed out that these periodic residuals can be modeled by a sum of the sinusoidal functions:
where n is the number of sinusoidal functions; and A_{k},f_{k}, and φ_{k} are the amplitude, frequency, and phase of the kth sinusoidal component, respectively.
3.2 Nonparametric regression
We propose an alternative residual correction method using a generalized additive model (GAM):
where s(∗) represents a smooth function. As an extension of the linear model by including additive smooth functions for the explanatory variables, the GAM has the potential to capture complex nonlinear patterns that a parametric model (e.g., highorder polynomials and sinusoidal functions) would miss.
The GAM software packages are widely available in various programming languages, such as R, Python, MATLAB, and SAS. Typical smooth functions include local polynomials, splines, Markov random fields, and Gaussian process smooths. In our experiments, the GAM regression of Eq. (21) was performed in R software using the “mgcv” package (Wood, 2022). A thinplate spline was chosen as the smoothing basis (i.e., the smooth function s), and the degree of smoothing was automatically determined by the generalized cross validation (GCV) criterion. For more on the theoretical foundations and technical details of the GAM method and the “mgcv” package, we refer the reader to Hastie and Tibshirani (1990) and Wood (2017).
4.1 Icesheet case study
4.1.1 Data processing
The comparative experiments of DEM coregistration and residual correction were carried out on 23 Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEM pairs from the western edge of the Greenland Ice Sheet (GrIS) (Fig. 3). Details of all ASTER DEM pairs are provided in the Supplement (Table S1), where two DEM pairs were used for visualization and analysis (Table 1). The raw stereoscopic DEMs were automatically produced by the US Geological Survey Land Processes Distributed Active Archive Center (LPDAAC) using SilcAst software (NASA et al., 2001).
The normalized difference bareness index (NDBI) was calculated from Landsat 8 images to extract stable regions (Nguyen et al., 2021):
where SWIR1 and GREEN represent the first shortwave infrared band (1.560–1.660 µm) and the green band (0.525–0.600 µm) of the Landsat 8 Operational Land Imager (OLI) data, respectively. All the terrainrelated information (slope, aspect, etc.), which served as explanatory variables of the regression, was then derived from the reference DEMs. In the coregistration and residual correction procedures, only DEM pixels over stable terrain were used for the regression, and a threesigma rule (i.e., more than 3 times the standard deviation) was employed on the elevation differences between two DEMs to remove erroneous data caused by misclassification of unstable terrain areas. A subset of the data (of no more than 50 000 pixels, to reduce the computational cost) was randomly selected as the training set, and the remaining pixels were used for the accuracy evaluation by comparing the median absolute difference (MedAD) (Mcmillan et al., 2019; Trevisani and Rocca, 2015):
where H_{Reference} and H_{Secondary} represent the reference and secondary DEM elevation, respectively.
4.1.2 DEM coregistration
Table 2 shows that all four coregistration methods effectively reduce the DEM biases, and the following findings were made by comparing the error statistics of the different algorithms.
The standard and linear versions of the NK method yield exactly the same outcomes. The only difference between the two algorithms is whether the regression equation is linear or not, which does not affect the coregistration results.
The NK simplified version produces similar results to the standard version. It should be noted that this conclusion may not hold true for other datasets, because it cannot be proven theoretically that approximating terrain slopes by their mean value would always lead to a reliable performance.
The RT method performs better than the three versions of the NK method. The coregistration errors of the RT method are smaller than those of the NK linear algorithm by an average of 4.6 % and a maximum of 15.3 %, which indicates that there are some scale and rotationinduced biases in the experimental DEM data.
Figure 4 shows the elevation differences of DEM pair GrIS1 before coregistration. All the pixels classified as water and potential outliers due to clouds were masked out for a better visualization, leaving the regions of bare land and glacier (bounded by the black lines). It can be seen from the figure that most pixels are negative values, indicating that the majority of the elevation differences are caused by vertical translation. Minor errors related to the terrain (induced by horizontal translation) and alongtrack coordinates (caused by jitter) can also be clearly observed.
The elevation difference maps (Fig. 5) demonstrate that the residuals of all three versions of the NK coregistration algorithms are consistent in terms of both magnitude and distribution. The RT algorithm shows better coregistration results, removing 11.8 % more errors compared to the NK linear version. A visual inspection of Fig. 5c and d reveals that the elevation differences of NK exhibit a positive trend in the northwest corner (the blue circle in Fig. 5c) and a negative trend in the southeast corner (the red circle in Fig. 5c), which are possibly caused by unconsidered attitude biases. In addition, some clustered outliers, which may consist of misclassified water and cloud pixels, can be clearly observed in the elevation difference maps (Fig. 5a). However, these outliers have little influence on the coregistration results because robust statistical methods (robust regression algorithms and a robust scale estimation method, i.e., MedAD) were used in the experiments.
4.1.3 Residual correction
The residual correction results for DEM pair GrIS1 are shown in Fig. 6. In the experiments, the polynomialfitting method used an eighthorder polynomial sequentially in the crosstrack and alongtrack directions, and the combination of the polynomial and the sumofsines method was implemented by first adopting an eighthorder polynomial in the crosstrack direction and then applying a sum of three sines in the alongtrack direction. A visual comparison reveals that the highorder polynomial removes the lowfrequency residuals only, whereas both the sum of sines and the GAM spline can capture the highfrequency signals. The MedAD values show that the GAM splinefitting method yields a higher accuracy than the two parametric regression methods.
The magnitude of the highfrequency signals in DEM pair GrIS2 is much greater than that in DEM pair GrIS1 in Fig. 7b. The polynomialfitting method again eliminates only the lowfrequency residuals. Figure 7c and d show that weak striped patterns exist in the residual results of both the sumofsines and the GAM splinefitting methods, indicating that the highfrequency errors are not completely removed. The MedAD values show that the combination of the polynomial and the sumofsines method is 5.1 % less accurate than the GAMbased method, which can be observed by the significant negative biases indicated by the arrows in Fig. 7c. Figure 8 further shows the fitting results in the alongtrack direction. The eighthorder polynomial only matches the longterm trend, and a followup experiment revealed that increasing the order of the polynomial still does not help to capture highfrequency signals. As marked by the red and purple arrows in Fig. 8 (corresponding to the regions indicated in Fig. 7), the maximum difference between the sumofsines and the GAM splinefitting results is about 5 m. Because the sinusoidal function is a parametric model whose parameters (amplitude, phase, and frequency) are global constants, there is no difference in shape between the different cycles. In contrast, the GAM spline yields a nonstrictly periodic curve by fitting the local relationship between the elevation differences (the response variable) and the alongtrack coordinates (the predictor variable) over parts of their range. A visual inspection shows that the GAM splinefitting results fit more closely with the local trends in the coregistration residuals, which indicates that the GAM splinefitting method might be a better alternative to the traditional parametric models for residual correction of DEM coregistration results, benefiting from its datadriven nature.
Finally, Table 3 summarizes the residual correction results for the 23 ASTER DEM pairs. The GAM splinefitting method outperforms the polynomial method and the combination of the polynomial and the sumofsines method by reducing 4.4 % and 2.1 % more residuals, respectively. We manually check the residual correction results for all the DEM pairs. A visual inspection shows that the remaining errors for a majority of the data (e.g., pair GrIS1 in Fig. 6) are almost randomly distributed in the scene. Only a few DEM pairs suffer from minor systematic errors caused by incompletely corrected jitter (e.g., pair GrIS2 in Fig. 7), where slight biases would be propagated into the glacier thickness change estimates.
4.2 Mountain glacier case study
4.2.1 Data processing
The mountain glacier experiments were performed on 22 DEM pairs from the Pamir region of High Mountain Asia (HMA) (Fig. 9), including ASTER DEMs, ZiYuan3 (ZY3) DEMs generated from ZiYuan3 tristereo optical scenes (Pan et al., 2013; Liu et al., 2020), and the global Shuttle Radar Topography Mission (SRTM) DEM (Farr et al., 2007) and the Copernicus DEM GLO30 (Copernicus, 2023) obtained using the Interferometric Synthetic Aperture Radar (InSAR) technique. Stable regions were extracted from three land cover classes (bare land, artificial surfaces, and cultivated land) in the GlobeLand30 land cover product (Jun et al., 2014; Li et al., 2021).
4.2.2 DEM coregistration
Like the icesheet case, the simplified and standard versions of the NK method yield similar results for the three test datasets of ZY3 DEMs, the SRTM DEM, and the Copernicus DEM in the HMA region (Table 5). The RT method shows better coregistration performance than the three versions of the NK method, removing 13.7 % more errors over the linear version on average.
Figure 10 depicts an example of the ZY3 DEM with large attitude errors. From Fig. 10a to c, it can be seen that the coregistration results of the NK method exhibit significant residuals in the southwest–northeast direction, leading to a false estimate of rapid glacier mass loss in the northern region. In contrast, the RT algorithm can effectively eliminate attitudeinduced bias and reduce coregistration errors by 83.3 % compared to the NK linear version. More test cases using ZY3 DEMs are provided in the Supplement, and their errorreduction ratios range from 0.7 % to 42.3 %. A visual comparison between Figs. 10d and 5d reveals that the coregistration residuals of the ZY3 DEM are much smaller than those of the ASTER DEM. This may be due to the fact that the ZY3 raw image has a resolution of 2.5–3.5 m (Zhang et al., 2018) and retains a high signaltonoise ratio after downsampling to 30 m.
Figure 11 illustrates the influence of the DEM coregistration method on glacier surface elevation change estimation. Panels a, b, and c show the coregistration results and glacier elevation change statistical values for an ASTER DEM pair acquired half a month apart (HMA2). Since the two DEMs were obtained at a very short time interval, it can be assumed that there is no elevation change. The result of the NK linear version contains obvious rotationinduced errors (Fig. 11a), the large proportion of missing data throughout the center of the image is caused by cloud cover, and the glaciercovered area in the southeast region only covers a small number of invalid pixels. We calculated the glacier elevation change values within the red circle and found a significant negative deviation (−6.396 m). No distinct abnormal trends were found in the results of the RT method, and the mean value of glacier elevation changes is very close to zero.
We provide DEM coregistration examples in the Supplement for more scenarios in HMA and New Zealand (NZL), such as a large number of glaciers, a large amount of vegetation, and a high noise level due to rough topography in the DEMs. Since no strong jitterinduced residuals were observed in the coregistration results of these DEM pairs, residual correction experiments were not performed.
The performance of different types of DEM coregistration methods has been intensively investigated by Paul et al. (2015) and Vacaflor et al. (2022). Their tests showed that the NK method achieved similar or better accuracy compared to many nonanalytical methods, such as the grid search method (Berthier et al., 2007), the LS3D method (Gruen and Akca, 2005), and the subwatershedbased method (Li et al., 2017), and therefore the NK method was recommended for practical applications due to the less computational effort (Paul et al., 2015).
This work focuses on the comparison of two analytical algorithms, the NK method and the RT method, which have been widely used in the cryosphere (Geyman et al., 2022; Hugonnet et al., 2021; Maurer et al., 2019) and photogrammetry (Aguilar et al., 2012; Zhang et al., 2010) studies, respectively. The characteristics of the two methods are summarized in Table 6, and the theoretical connections and differences between them are discussed in the following.

The form of the regression.
The NK method can be expressed as either a nonlinear or linear equation, while the RT method only employs a linear regression model. The disadvantages of nonlinear regression over linear regression are that it works iteratively and often requires more computational resources.

The explanatory variables in the regression.
The NK method was inspired by the similarity between an elevation difference map and a hillshade, which is predicted based on the terrain slope and aspect. The RT method, on the other hand, employs the terrain gradients (i.e., the partial first derivatives) in the X and Y directions as the explanatory variables. From Eqs. (15) and (16), it can be seen that the two groups of terrain variables are actually equivalent.

Regression coefficients.
In the RT method, the misalignment between two DEMs is modeled by a 3D similarity transformation, including three translation, one scale, and three rotation factors. The NK method considers the spatial shift only, and the regression coefficients can be either cylindrical coordinates ($a,b,c$) or Cartesian coordinates ($\mathrm{\Delta}X,\mathrm{\Delta}Y,\mathrm{\Delta}Z$) of the shift vector.
Based on the above analysis, it can be concluded that the number of regression coefficients is the only significant difference between the NK and RT methods. In other words, the RT method can be viewed as an extension of NK by additionally modeling the scale and rotation errors. The advantage of the original NK method (i.e., the simplified version) is that only one explanatory variable (ψ) exists in the regression model, and the shift vector can therefore be calculated by a curvefitting technique or estimated from a scatter plot, which is easy to adopt for users with a limited knowledge of statistics. However, in terms of precision, the NK method is theoretically inferior to the RT method, because only the shiftinduced errors are considered.
The NH method (Noh and Howat, 2014) is another one analytical algorithm that has been previously used in glacial studies, in addition to the NK method. It can be theoretically proven that the NH method is basically equivalent to the RT method (details are available in the Supplement), because both of them are derived from a 3D similarity transformation considering the errors induced by three translation, one scale, and three rotation factors. The disadvantage of the NH method is that small Euler angles are not approximated, and, accordingly, the regression equation has a very complicated form (see Eqs. S3 and S4 in the Supplement), which hinders the replication of the algorithm by other researchers.
To sum up, the three analytical algorithms for the DEM coregistration problem, i.e., NK, NH, and RT, have a strong theoretical relationship, despite being presented in diverse forms in the original literature. Taking into account both the algorithm accuracy and ease of use in practical applications, we recommend applying the RT method instead of the NK and NH methods in glacial studies.
A residual correction procedure is required after DEM coregistration if satellite attitude jitters and other complex errors are not be properly removed before the production of DEMs. The residual signals often appear at several frequencies. The numbers are not constant for different DEM sources (Girod et al., 2017; Wang et al., 2016; Ye et al., 2019), and they primarily depend upon sensors and data preprocessing techniques. For the traditional parametric regression methods, the optimization choice of the degree of polynomials (or the number of sinusoidal functions) takes n trials, where n is the highest degree of the polynomial (or the maximum number of sines) allowed. For example, n=6 was suggested by Girod et al. (2017). If the residual signals in the data do not have a repeated regular shape, the performance of traditional parametric regression methods is limited by their predefined models. The GAM splinefitting method is datadriven and can capture more complex nonlinear patterns in residuals. The spline parameters can be automatically determined by the generalized cross validation or other criterion, and the optimization process requires high computational cost. The GAM spline fitting can be used as an alternative to traditional parametric regression models when the residuals cannot be precisely fitted by highorder polynomials or a sum of sines. It should be noted that if the noise level of DEM coregistration residuals is very high, the performance of the GAM splinefitting method is just slightly better than that of traditional parametric regression methods.
In Sect. 4, we focused on the accuracy comparison of various algorithms and conducted tests on DEM pairs with good geometric conditions. However, the geometric constraint of stable terrain may be weak for the DEM pairs located at the edge of an ice sheet with little stable terrain or covered by heavy clouds. Under such conditions, the reliability issue of algorithms becomes even more critical.
Figure 12 shows a representative example of one ASTER DEM pair located on the western edge of GrIS, where the stable terrain is geographically distributed in the southwest corner only. The time interval between reference DEM (ASTER DEM 20190725) and secondary DEM (ASTER DEM 20190826) is 1 month, and therefore the ice surface elevation can be considered unchanged. Although the RT algorithm yields significantly smaller coregistration residuals in the stable region than the NK method, it is prone to producing larger biases over icecovered regions. These biases cannot be removed by the residual correction procedure, because the residual trend over icecovered regions is completely different from that over stable regions.
The second example is the residual regression of the DEM pair GrIS19 (in Table S1) by taking the terrain elevation as the explanatory variable. It can be seen from Fig. 13 that the mean elevation of glaciers is much higher than that of bare lands, and a long extrapolation is therefore required. The prediction results of the polynomial and splinefitting methods are strongly biased in highaltitude regions (>500 m).
It can be noticed from the above two examples that, in the DEM coregistration and residual correction tasks, unreliable regression results cannot be detected from the elevation differences of stable regions. As a rule of thumb, when a data extrapolation is needed, it is recommended to adopt a conservative strategy by decreasing the degree of freedom of the regression model, e.g., dropping some explanatory variables (in DEM coregistration), and reducing the degree of the polynomial (in highorder polynomial regression) or smoothing (in spline regression).
Many researchers pointed out the need for residual correction along terrain altitudes after DEM coregistration. This correction is absolutely necessary for InSAR DEMs to reduce the biases caused by radar penetration (Gardelle et al., 2012; Li et al., 2021), while it is only optional for stereoscopic DEMs (Nuth and Kääb, 2011) because there is no strong physical basis to explain it. Given that obvious elevationdependent biases were not observed in our experiments (e.g., Fig. 13a), the terrain elevation was not introduced as an explanatory variable in our residual regression model (i.e., its degree of freedom is zero), and, accordingly, the problem of data extrapolation along terrain altitudes was circumvented.
In this paper, we have made a thorough comparison of the DEM coregistration methods of Nuth and Kääb (2011) and Rosenholm and Torlegard (1988), and we proposed a GAMbased method to correct DEM coregistration residuals. The theoretical analysis and experimental results support the following conclusions.

The NK method and the RT method are theoretically compatible with each other. On the one hand, the terrainrelated information used by the two methods as explanatory variables in their regressions – slope/aspect and gradient – can be proven to be equivalent through theoretical analysis. On the other hand, even though the methods of NK and RT utilize distinct regression forms, the nonlinear regression equation used by the former can be converted into a linear equation with the same structure as the latter.

Rotation and scale biases should be taken into account in DEM coregistration. The only significant difference between the methods of NK and RT is that the latter models the translation, scale, and rotationinduced biases, while the former only considers the spatial translation. Comparative experiments conducted on multiple DEM pairs showed that the RT method consistently outperformed the NK method in terms of coregistration residuals. Thus, we recommend applying the RT method also in glacial studies.

GAM spline fitting can be used as an alternative to traditional parametric regression models in correcting DEM coregistration residuals. ASTER DEMs often suffer from some complex errors with multiple frequencies induced by satellite attitude jitter. Benefiting from its datadriven nature, the GAM splinefitting method can capture the complex nonlinear patterns in DEM coregistration residuals, whereas the performance of the parametric regression methods is limited by their predefined models when the residuals do not have a repeated regular shape.
The MATLAB code for DEM coregistration is available at https://github.com/shenapm/DemCoReg (last access: 22 October 2023) and archived at https://doi.org/10.5281/zenodo.8098337 (Li and Shen, 2023).
ASTER DEMs are freely available at https://doi.org/10.5067/ASTER/AST14DEM.003 (NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team, 2001). Landsat 8 images used in this study can be obtained at no cost from the United States Geological Survey (https://earthexplorer.usgs.gov/, last access: 6 December 2023). The Greenland drainage basins can be obtained from http://imbie.org/imbie3/drainagebasins/ (last access: 6 December 2023; Rignot and Mouginot, 2012; IMBIE, 2023).
The supplement related to this article is available online at: https://doi.org/10.5194/tc1752992023supplement.
XS and TL conceptualized and initiated the study. TL performed the data processing and analyses, prepared the figures and tables, and wrote the draft manuscript. XS contributed to the review and improvement of the manuscript. YH, BL, LJ, and HW assisted in the editing and refining of the manuscript.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We are grateful to Tobias Bolch, Anton Schenk, and an anonymous referee, whose insight and suggestions improved the paper.
This research was supported by the National Key R&D Program of China (nos. 2017YFA0603103 and 2018YFC1406102); the Strategic Priority Research Program of the Chinese Academy of Sciences (no. XDA19070104); the National Natural Science Foundation of China (nos. 41771490, 41974009, 42004007 and 42174046); the State Key Laboratory of Geodesy and Earth's Dynamics, Innovation Academy for Precision Measurement Science and Technology, CAS (nos. S21L6102 and S22L6101); and the Key Research Program of Frontier Sciences, CAS (nos. QYZDBSSWDQC042 and QYZDJSSWDQC027).
This paper was edited by Tobias Bolch and reviewed by Anton Schenk and one anonymous referee.
Aguilar, F. J., Aguilar, M. A., Fernandez, I., Negreiros, J. G., Delgado, J., and Perez, J. L.: A New TwoStep Robust Surface Matching Approach for ThreeDimensional Georeferencing of Historical Digital Elevation Models, IEEE Geosci. Remote S., 9, 589–593, https://doi.org/10.1109/LGRS.2011.2175899, 2012.
Akca, D.: Coregistration of Surfaces by 3D Least Squares Matching, Photogramm. Eng. Rem. S., 76, 307–318, https://doi.org/10.14358/PERS.76.3.307, 2010.
Berthier, E., Cabot, V., Vincent, C., and Six, D.: Decadal RegionWide and GlacierWide Mass Balances Derived from MultiTemporal ASTER Satellite Digital Elevation Models. Validation over the MontBlanc Area, Front. Earth Sci., 4, 63, https://doi.org/10.3389/feart.2016.00063, 2016.
Berthier, E., Arnaud, Y., Kumar, R., Ahmad, S., Wagnon, P., and Chevallier, P.: Remote sensing estimates of glacier mass balances in the Himachal Pradesh (Western Himalaya, India), Remote Sens. Environ., 108, 327–338, https://doi.org/10.1016/j.rse.2006.11.017, 2007.
Besl, P. J. and Mckay, N. D.: A Method for Registration of 3D Shapes, IEEE T. Pattern Anal., 14, 239–256, https://doi.org/10.1109/34.121791, 1992.
Bolch, T., Pieczonka, T., and Benn, D. I.: Multidecadal mass loss of glaciers in the Everest area (Nepal Himalaya) derived from stereo imagery, The Cryosphere, 5, 349–358, https://doi.org/10.5194/tc53492011, 2011.
Brun, F., Berthier, E., Wagnon, P., Kääb, A., and Treichler, D.: A spatially resolved estimate of High Mountain Asia glacier mass balances from 2000 to 2016, Nat. Geosci., 10, 668–673, https://doi.org/10.1038/ngeo2999, 2017.
Copernicus: Copernicus DEM – Global and European Digital Elevation Model (COPDEM), GLO30, ESA, Copernicus [data set], https://doi.org/10.5270/ESAc5d3d65, 2023.
Cucchiaro, S., Maset, E., Cavalli, M., Crema, S., Marchi, L., Beinat, A., and Cazorzi, F.: How does coregistration affect geomorphic change estimates in multitemporal surveys?, GISci. Remote Sens., 57, 611–632, https://doi.org/10.1080/15481603.2020.1763048, 2020.
Di, K., Hu, W., Liu, Y., and Peng, M.: Coregistration of Chang'E1 stereo images and laser altimeter data with crossover adjustment and image sensor model refinement, Adv. Space Res., 50, 1615–1628, https://doi.org/10.1016/j.asr.2012.06.037, 2012.
Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., Kobrick, M., Paller, M., Rodriguez, E., Roth, L., Seal, D., Shaffer, S., Shimada, J., Umland, J., Werner, M., Oskin, M., Burbank, D., and Alsdorf, D.: The Shuttle Radar Topography Mission, Rev. Geophys., 45, RG2004, https://doi.org/10.1029/2005RG000183, 2007.
Gardelle, J., Berthier, E., and Arnaud, Y.: Impact of resolution and radar penetration on glacier elevation changes computed from DEM differencing, J. Glaciol., 58, 419–422, https://doi.org/10.3189/2012JoG11J175, 2012.
Gardelle, J., Berthier, E., Arnaud, Y., and Kääb, A.: Regionwide glacier mass balances over the PamirKarakoramHimalaya during 1999–2011, The Cryosphere, 7, 1263–1286, https://doi.org/10.5194/tc712632013, 2013.
Geyman, E. C., van Pelt, W. J. J., Maloof, A. C., Aas, H. F., and Kohler, J.: Historical glacier change on Svalbard predicts doubling of mass loss by 2100, Nature, 601, 374–379, https://doi.org/10.1038/s41586021043144, 2022.
Girod, L., Nuth, C., Kääb, A., McNabb, R., and Galland, O.: MMASTER: Improved ASTER DEMs for Elevation Change Monitoring, Remote Sens., 9, 704, https://doi.org/10.3390/rs9070704, 2017.
Gorokhovich, Y. and Voustianiouk, A.: Accuracy assessment of the processed SRTMbased elevation data by CGIAR using field data from USA and Thailand and its relation to the terrain characteristics, Remote Sens. Environ., 104, 409–415, https://doi.org/10.1016/j.rse.2006.05.012, 2006.
Gruen, A. and Akca, D.: Least squares 3D surface and curve matching, ISPRS J. Photogramm., 59, 151–174, https://doi.org/10.1016/j.isprsjprs.2005.02.006, 2005.
Hastie, T. and Tibshirani, R.: Generalized additive models, Chapman and Hall, London, https://doi.org/10.1201/9780203753781, 1990.
Hofton, M., Dubayah, R., Blair, J. B., and Rabine, D.: Validation of SRTM elevations over vegetated and nonvegetated terrain using medium footprint lidar, Photogramm. Eng. Rem. S., 72, 279–285, https://doi.org/10.14358/PERS.72.3.279, 2006.
Hugonnet, R., McNabb, R., Berthier, E., Menounos, B., Nuth, C., Girod, L., Farinotti, D., Huss, M., Dussaillant, I., Brun, F., and Kääb, A.: Accelerated global glacier mass loss in the early twentyfirst century, Nature, 592, 726–731, https://doi.org/10.1038/s4158602103436z, 2021.
IMBIE: IMBIE3 Rignot Greenland Drainage Basins, IMBIE [data set], http://imbie.org/imbie3/drainagebasins/, last access: 6 December 2023.
Jun, C., Ban, Y., and Li, S.: Open access to Earth landcover map, Nature, 514, 434–434, https://doi.org/10.1038/514434c, 2014.
Karkee, M., Steward, B. L., and Aziz, S. A.: Improving quality of public domain digital elevation models through data fusion, Biosyst. Eng., 101, 293–305, https://doi.org/10.1016/j.biosystemseng.2008.09.010, 2008.
Kim, T. and Jeong, J.: DEM matching for bias compensation of rigorous pushbroom sensor models, ISPRS J. Photogramm., 66, 692–699, https://doi.org/10.1016/j.isprsjprs.2011.06.002, 2011.
Li, C., Jiang, L., Liu, L., and Wang, H.: Regional and AltitudeDependent Estimate of the SRTM C/XBand Radar Penetration Difference on High Mountain Asia Glaciers, IEEE J. Sel. Top. Appl., 14, 4244–4253, https://doi.org/10.1109/jstars.2021.3070362, 2021.
Li, H., Deng, Q. L., and Wang, L. C.: Automatic CoRegistration of Digital Elevation Models Based on Centroids of Subwatersheds, IEEE T. Geosci. Remote, 55, 6639–6650, https://doi.org/10.1109/Tgrs.2017.2731048, 2017.
Li, T. and Shen, X.: Coregistration and residual correction of digital elevation models: A comparative study, Zenodo [code], https://doi.org/10.5281/zenodo.8098337, 2023.
Liu, L., Jiang, L., Jiang, H., Wang, H., Ma, N., and Xu, H.: Accelerated glacier mass loss (2011–2016) over the Puruogangri ice field in the inner Tibetan Plateau revealed by bistatic InSAR measurements, Remote Sens. Environ., 231, 111241, https://doi.org/10.1016/j.rse.2019.111241, 2019.
Liu, L., Jiang, L., Zhang, Z., Wang, H., and Ding, X.: Recent Accelerating Glacier Mass Loss of the Geladandong Mountain, Inner Tibetan Plateau, Estimated from ZiYuan3 and TanDEMX Measurements, Remote Sens., 12, 472, https://doi.org/10.3390/rs12030472, 2020.
Maurer, J. M., Schaefer, J. M., Rupper, S., and Corley, A.: Acceleration of ice loss across the Himalayas over the past 40 years, Sci. Adv., 5, eaav7266, https://doi.org/10.1126/sciadv.aav7266, 2019.
McMillan, M., Muir, A., Shepherd, A., Escolà, R., Roca, M., Aublanc, J., Thibaut, P., Restano, M., Ambrozio, A., and Benveniste, J.: Sentinel3 DelayDoppler altimetry over Antarctica, The Cryosphere, 13, 709–722, https://doi.org/10.5194/tc137092019, 2019.
NASA, METI, AIST, Japan Spacesystems and U.S./Japan ASTER Science Team: ASTER DEM Product, NASA EOSDIS Land Processes DAAC [data set], https://doi.org/10.5067/ASTER/AST14DEM.003, 2001.
Nguyen, C. T., Chidthaisong, A., Kieu Diem, P., and Huo, L.Z.: A Modified Bare Soil Index to Identify Bare Land Features during Agricultural FallowPeriod in Southeast Asia Using Landsat 8, Land, 10, 231, https://doi.org/10.3390/land10030231, 2021.
Noh, M. and Howat, I. M.: Automated Coregistration of Repeat Digital Elevation Models for Surface Elevation Change Measurement Using Geometric Constraints, IEEE T. Geosci. Remote, 52, 2247–2260, https://doi.org/10.1109/TGRS.2013.2258928, 2014.
Nuth, C. and Kääb, A.: Coregistration and bias corrections of satellite elevation data sets for quantifying glacier thickness change, The Cryosphere, 5, 271–290, https://doi.org/10.5194/tc52712011, 2011.
Pan, H., Zhang, G., Tang, X., Li, D., Zhu, X., Zhou, P., and Jiang, Y.: Basic Products of the ZiYuan3 Satellite and Accuracy Evaluation, Photogramm. Eng. Rem. S., 79, 1131–1145, https://doi.org/10.14358/PERS.79.12.1131, 2013.
Paul, F., Bolch, T., Kääb, A., Nagler, T., Nuth, C., Scharrer, K., Shepherd, A., Strozzi, T., Ticconi, F., Bhambri, R., Berthier, E., Bevan, S., Gourmelen, N., Heid, T., Jeong, S., Kunz, M., Lauknes, T. R., Luckman, A., Merryman Boncori, J. P., Moholdt, G., Muir, A., Neelmeijer, J., Rankl, M., VanLooy, J., and Van Niel, T.: The glaciers climate change initiative: Methods for creating glacier area, elevation change and velocity products, Remote Sens. Environ., 162, 408–426, https://doi.org/10.1016/j.rse.2013.07.043, 2015.
Peckham, R. J. and Jordan, G.: Digital Terrain Modelling: Development and Applications in a Policy Support Environment, Springer Berlin, Heidelberg, https://doi.org/10.1007/9783540367314, 2007.
Peduzzi, P., Herold, C., and Silverio, W.: Assessing high altitude glacier thickness, volume and area changes using field, GIS and remote sensing techniques: the case of Nevado Coropuna (Peru), The Cryosphere, 4, 313–323, https://doi.org/10.5194/tc43132010, 2010.
Pfeffer, W. T., Arendt, A. A., Bliss, A., Bolch, T., Cogley, J. G., Gardner, A. S., Hagen, J.O., Hock, R., Kaser, G., Kienholz, C., Miles, E. S., Moholdt, G., Mölg, N., Paul, F., Radić, V., Rastner, P., Raup, B. H., Rich, J., Sharp, M. J., and The Randolph Consortium: The Randolph Glacier Inventory: a globally complete inventory of glaciers, J. Glaciol., 60, 537–552, https://doi.org/10.3189/2014JoG13J176, 2014.
Pieczonka, T., Bolch, T., Wei, J., and Liu, S.: Heterogeneous mass loss of glaciers in the AksuTarim Catchment (Central Tien Shan) revealed by 1976 KH9 Hexagon and 2009 SPOT5 stereo imagery, Remote Sens. Environ., 130, 233–244, https://doi.org/10.1016/j.rse.2012.11.020, 2013.
Rastner, P., Bolch, T., Mölg, N., Machguth, H., Le Bris, R., and Paul, F.: The first complete inventory of the local glaciers and ice caps on Greenland, The Cryosphere, 6, 1483–1495, https://doi.org/10.5194/tc614832012, 2012.
RGI Consortium: Randolph Glacier Inventory – A Dataset of Global Glacier Outlines, Version 6, Boulder, Colorado USA, National Snow and Ice Data Center [data set], https://doi.org/10.7265/4m1fgd79, 2017.
Rignot, E. and Mouginot, J.: Ice flow in Greenland for the International Polar Year 2008–2009, Geophys. Res. Lett., 39, L11501, https://doi.org/10.1029/2012GL051634, 2012.
Rodriguez, E., Morris, C. S., and Belz, J. E.: A global assessment of the SRTM performance, Photogramm. Eng. Rem. S., 72, 249–260, https://doi.org/10.14358/Pers.72.3.249, 2006.
Rosenholm, D. and Torlegard, K.: Threedimensional absolute orientation of stereo models using digital elevation models, Photogramm. Eng. Rem. S., 54, 1385–1389, 1988.
Rusinkiewicz, S. and Levoy, M.: Efficient variants of the ICP algorithm, Proceedings third international conference on 3D digital imaging and modeling, Quebec City, Canada, 28 May–1 June 2001, 145–152, https://doi.org/10.1109/IM.2001.924423, 2001.
Sedaghat, A. and Naeini, A. A.: DEM orientation based on local feature correspondence with global DEMs, GISci. Remote Sens., 55, 110–129, https://doi.org/10.1080/15481603.2017.1364879, 2018.
Trevisani, S. and Rocca, M.: MAD: robust image texture analysis for applications in high resolution geomorphometry, Comput. Geosci., 81, 78–92, https://doi.org/10.1016/j.cageo.2015.04.003, 2015.
United States Geological Survey: USGS EarthExplorer [data set], https://earthexplorer.usgs.gov/, last access: 6 December 2023.
Vacaflor, P., Lenzano, M. G., Vich, A., and Lenzano, L.: CoRegistration Methods and Error Analysis for Four Decades (1979–2018) of Glacier Elevation Changes in the Southern Patagonian Icefield, Remote Sens., 14, 820, https://doi.org/10.3390/rs14040820, 2022.
Van Niel, T. G., McVicar, T. R., Li, L. T., Gallant, J. C., and Yang, Q. K.: The impact of misregistration on SRTM and DEM image differences, Remote Sens. Environ., 112, 2430–2442, https://doi.org/10.1016/j.rse.2007.11.003, 2008.
Wang, M., Zhu, Y., Jin, S. Y., Pan, J., and Zhu, Q. S.: Correction of ZY3 image distortion caused by satellite jitter via virtual steady reimaging using attitude data, ISPRS J. Photogramm., 119, 108–123, https://doi.org/10.1016/j.isprsjprs.2016.05.012, 2016.
Wood, S. N.: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation, R package “mgcv”, version 1.840 [code], https://cran.rproject.org/web/packages/mgcv/index.html (last access: 6 December 2023), 2022.
Wood, S. N.: Generalized Additive Models: An Introduction with R, Second Edition, Texts in Statistical Science, Chapman & Hall/CRC, https://doi.org/10.1201/9781315370279, 2017.
Wu, B., Guo, J., Hu, H., Li, Z. L., and Chen, Y. Q.: Coregistration of lunar topographic models derived from Chang'E1, SELENE, and LRO laser altimeter data based on a novel surface matching method, Earth Planet. Sc. Lett., 364, 68–84, https://doi.org/10.1016/j.epsl.2012.12.024, 2013.
Ye, Z., Xu, Y. S., Tong, X. H., Zheng, S. Z., Zhang, H., Xie, H., and Stilla, U.: Estimation and analysis of alongtrack attitude jitter of ZiYuan3 satellite based on relative residuals of triband multispectral imagery, ISPRS J. Photogramm., 158, 188–200, https://doi.org/10.1016/j.isprsjprs.2019.10.012, 2019.
Zhang, T., Cen, M., Ren, Z., Yang, R., Feng, Y., and Zhu, J.: Ability to detect and locate gross errors on DEM matching algorithm, Int. J. Digit. Earth, 3, 72–82, https://doi.org/10.1080/17538940903033175, 2010.
Zhang, T., Lei, B., Wang, J., Li, Y., Liu, K., and Li, T.: Preliminary Quality Analysis of the Triple LinearArray and Mul Tispectral Images of ZY3 02 Satellite, in: Proceedings of the IEEE International Geoscience and Remote Sensing Symposium, Valencia, Spain, 23–27 July 2018, 9172–9175, https://doi.org/10.1109/IGARSS.2018.8518600, 2018.