ABSTRACT
ABSTRACT
Digital elevation models (DEMs) represent the Earth’s topography and support a variety of applications, ranging from extracting watershed drainage structure to measuring glacier volume. Applications that require conflation of multiple DEMs are problematic, however, because of misregistration. We explore a spatial optimization technique to quantify the misregistration on the pixel level between NASA’s Shuttle Radar Topography Mission (SRTM) elevation data set and the USGS’s National Elevation Dataset (NED). The misregistration, estimated within blocks by horizontal offset and direction, is modelled spatially in terms of typical topographical parameters: slope, aspect and elevation. A Nelder–Mead algorithm was implemented to compute the local shift at two study sites in the San Gabriel Mountains (SGM) and Santa Monica Mountains (SMM) in Los Angeles County, California. The magnitude of misregistration is generally less than the distance between DEM postings, and varies by terrain type, being larger in steeper terrains; nevertheless, misregistration is systematic within these continuous terrains that have similar topographical features.
1. Introduction
Recent years have seen a proliferation of digital elevation data, largely as a result of new acquisition systems (Thompson, Bell, and Butler Citation2001). Digitized contours are often obtained by scanning published topographic maps, or acquired as an intermediate product in the mapping process. Rasters of elevation data are often used in geographic applications and analysis. These data sets are often found in the form of digital elevation models (DEMs), through interpolation from contours (the contour-to-grid or CTOG process), by remote sensing using radar (e.g., Shuttle Radar Topography Mission or SRTM data) or LiDAR data, or by interpolation from spot heights. The resulting data sets will vary by spatial resolution and the accuracy of position and elevation, and also by semantics, or precisely what is meant by height. LiDAR first returns, for example, will record the elevation of the tree canopy or buildings; radar will also be affected by dense canopies and buildings (DEMs acquired from first-return LiDAR or radar are sometimes termed digital surface models or DSMs for this reason); and the elevations depicted on topographic maps are typically those of the bare ground devoid of buildings, trees, and in some cases sand dunes (Zhang and Whitman Citation2005) (bare-earth DEMs obtained from topographic maps are sometimes termed digital terrain models, or DTMs). Some DEMs record the elevation measured at point locations, while in other cases the elevation reported may be the average over a pixel, or the result of a more complex spatial convolution.
This variety of elevation data is opening new applications that rely on the ability to fuse or compare two or more sources. Our study, for example, is motivated by a need to estimate the spatial distribution of total building volume in Los Angeles at a coarse scale, and our exploration of the possibility to accomplish this by comparing a DEM that records the elevation of bare ground with one, specifically SRTM that records the top surface of any buildings that are present. This effort required a careful registration of the two DEMs; this article reports on the technique we developed to address that problem. The broader issue of building volume estimation will be addressed in a separate article. Besides the application of DEMs in our study, a number of other scientific research topics, such as the extraction of watershed drainage structure (Maddalena Citation2010), estimation of canopy height (Miliaresis and Delikaraoglou Citation2009), and measurement of glacier volume (Nuth and Kääb Citation2011), all need cross-source DEM analysis. To perform precise science, it is imperative that the misregistration across DEMs is identified and as far as possible corrected.
In some cases, it is possible to address misregistration visually, or at least to gain some rough impression of its amount and direction. For example, in the Ridge and Valley Province of the Appalachian Mountains, an offset between two DEMs that is perpendicular to the ridges would lead to a characteristic pattern that alternates positive and negative elevation errors. Figure 1(a) shows an area near Harrisonburg, Virginia, where the ridges run approximately southwest to northeast. We then simulated a misregistration of 0.83 pixel perpendicular to the ridges. Figure 1(b) shows the result. Along each mountain ridge, the differences are all negative (black) on the northwest facing slope of the ridge, and positive on the southeast slope, allowing the ridges and valleys to become clear, as if the surface were illuminated from the southeast.
Published online:
10 September 2015Figure 1. (a) Satellite imagery of Appalachian Mountain ranges near Harrisonburg, Virginia. (b) Simulation of the effects of a 0.83-pixel misregistration.
Figure 1. (a) Satellite imagery of Appalachian Mountain ranges near Harrisonburg, Virginia. (b) Simulation of the effects of a 0.83-pixel misregistration.
The Appalachian Mountains example is an ideal case illustrating the patterns caused by subpixel misregistration because (1) these mountains are mostly single-directional; (2) we impose a constant shift in a single direction perpendicular to the structure of the mountain range. Therefore, the misregistration can be detected easily. However, in real-world situations, the identification of misregistration is much more complicated. It is therefore essential to develop a computer-aided spatial model for automated misregistration detection. In this article, we report on our efforts to design a spatial model that considers three important topographic parameters: slope, aspect and elevation, in order to quantify the misregistration (we also term it the shift) in hilly areas using optimization techniques.
The rest of the paper is organized as follows. Section 2 introduces the related research in the literature; Sections 3–5 describe the study area, the data set in use and the proposed method; Section 6 presents experimental results from our case studies and Section 7 concludes the article and identifies future research directions.
2. Literature
In image processing, it is sometimes necessary to fuse two or more images that have been produced by different approaches, acquired at different times or acquired from different points of view (Fonseca and Manjunath Citation1996; Zitová and Flusser Citation2003). In such cases, it is desirable to remove any misregistration that may be present, using one of two traditional methods. The feature-matching method involves computing the cross-correlation for a range of distances and directions between the image to be registered and a reference image, and selecting the distance and direction for which cross-correlation is maximum (Reddy and Chatterji Citation1996; Habib et al. Citation2005; Habib and Al-Ruzouq Citation2005; Guizar-Sicairos, Thurman, and Fienup Citation2008). Transform-model estimation involves the identification of a number of control points on both the image to be registered and the reference image that can be used to warp one to fit the other. Global approaches such as shape-preserving mapping and local approaches such as thin-plate splines (Davis et al. Citation1997) are commonly used transform models.
Most of the image-registration studies reported in the literature focused on registering images that exhibit sharp changes of value. Land-cover images, for example, have values that change sharply at the edges of different land-cover types, and may be classified, in which case values are measured on nominal scales. In contrast, elevation surfaces are measured on interval scales and typically follow Tobler’s first law of geography (Sui Citation2004), exhibiting strong spatial autocorrelations over short distances.
One popular method of measuring the positional accuracy of a DEM and improving its registration is by using ground control points (GCPs). The basic procedure is to: (1) identify a set of GCPs and measure their locations; (2) obtain the elevation at each GCP from in-situ measurement or from an independent source that is known to have better accuracy (Maune, Maitra, and McKay Citation2007); (3) identify feature points to be matched on the DEM and (4) estimate the parameters of a 3D transformation using the matched point pairs, with the objective of minimizing the least-square error of the elevations. However, this method is not effective for large-scale terrain DEM registration because it is costly and time consuming, and it may occasionally be impossible to obtain enough GCPs (Goodchild, Buttenfield, and Wood Citation1994; Zebker et al. Citation1994). Moreover, if the number of GCPs is limited and the transform model is poorly chosen, the result may be large discrepancies between the transformed DEM and the reference source, even after transformation (Carpenter and Hogarty, Citation2007).
To overcome this drawback, researchers have investigated the capabilities of DEM-to-DEM registration utilizing all of the available terrain information (Carpenter and Hogarty, Citation2007; Dawn, Saxena, and Sharma Citation2010). Carpenter and Hogarty (Citation2007) performed a simultaneous vertical and horizontal DEM registration using a ‘coarse exhaustive’ algorithm, a ‘coarse with parabolic refinement’ algorithm and a ‘unified least squares’ algorithm. The objective was to find the best fit between the master DEM (finer-resolution DEM) and the DEM to be registered (coarser-resolution DEM). This study found that the shift obtained is more reliable using DEMs with similar resolutions than with very different resolutions. It was not mentioned whether the complexity of terrains, such as local slope and aspect, were considered in the registration.
Recently, subpixel misregistration of images has attracted much attention. Dai and Khorram (Citation1998) showed that highly accurate change detection (with less than 10% error induced) based on multi-temporal Landsat Thematic Mapper (TM) images requires that the magnitude of misregistration be less than 0.2 pixel (see also Townshend et al. Citation1992), in other words, less than 0.2 times the horizontal separation of DEM postings. Simard et al. (Citation2006) reported an instance of subpixel DEM misregistration. Van Niel et al. (Citation2008) investigated the impact of misregistration between DEMs, including SRTM, and found that differences between DEMs are very sensitive to subpixel misregistration, and that the sensitivity depends on the pixel sizes of the data sets being compared (see also Smith and Sandwell Citation2003). They also concluded that even low levels of misregistration caused considerable differences in landscapes with steep terrains. However, neither of the above studies quantified the spatial variation of misregistration and its relationship to terrain characteristics. To extend previous studies, Streutker, Glenn, and Shrestha (Citation2011) proposed a method to utilize slope for the co-registration of overlapping elevation surface. Nuth and Kääb (Citation2011) modelled the misregistration among SRTM, ICESat and the ASTER GDEM data sets by considering elevation, local slope and aspect. This is a pioneer work in utilizing a complete set of topographical parameters to quantify the misregistration, however, we argue that horizontal shift has local characteristic instead of having a globally consistent offset. This spatial variation is significant especially in hilly regions where the changes of terrains vary by subareas. Additionally, as the goal of Nuth and Kääb’s research is to measure the thickness change of glacier, substantial implementation details on the co-registration are missing. Our method, in comparison, extends the prior work of Nuth and Kääb (Citation2011), and obtains block-by-block estimates of misregistration between two DEMs, allowing us to analyse the relationship between misregistration, location and the characteristics of terrain using a fine-tuned optimization approach. Once the misregistration is corrected, any systematic differences in elevations between the two data sets can be evaluated (Rodríguez, Morris, and Belz Citation2006).
3. Study sites
Two study sites in Los Angeles County in Southern California area were selected for this study: one located in the San Gabriel Mountains in the north of the county, and the other in the Santa Monica Mountains in the southwest, as Figure 2(a) shows. These two mountainous regions were selected because they contain few buildings and tree cover is sparse, and therefore the elevation difference between DSMs, such as SRTM, and DTMs, such as NED (National Elevation Dataset), would be minimal. The study site SGM (black box in Figure 2(b)) in the San Gabriel Mountain region covers 153 km2 [−118.38°, 34.34°; −118.19°, 34.42°]. It has very steep terrain, the distribution of slopes [slope defined as vertical rise over horizontal run, that is, the tangent of the slope angle, and estimated at DEM points using the Horn (Citation1981) algorithm] having a mean of 50% and a standard deviation of 22%. The elevations in this area lie between 568 m and 1845 m (mean = 1145 m; standard deviation = 249.5 m). The study site SMM (black box in Figure 2(c)) in the Santa Monica Mountain region covers 274 km2 [−118.86°, 34.04°;-118.56°, 34.13°]. Although less steep than SGM, the SMM region also has relatively steep terrain, with a mean slope of 38% (standard deviation 21%). The maximum and minimum elevations of this site are 3.2 m and 860 m, respectively. The mean elevation is 374 m, with a standard deviation of 146.8 m.
Published online:
10 September 2015Figure 2. (a) The regional placement of the study sites. (b, c) The topography of the San Gabriel Mountain region (b) and the Santa Monica Mountain region (c) from the NED DEM (the black rectangles indicate the study areas).
Figure 2. (a) The regional placement of the study sites. (b, c) The topography of the San Gabriel Mountain region (b) and the Santa Monica Mountain region (c) from the NED DEM (the black rectangles indicate the study areas).
4. Data
The Shuttle Radar Topography Mission (Farr et al. Citation2007), an international effort led by the U.S. National Geospatial-Intelligence Agency and the U.S. National Aeronautics and Space Administration, was carried out between 11 February and 22 February 2000. A C/X-band Synthetic Aperture Radar (SAR) was carried on board the Space Shuttle Endeavour, and used to scan the Earth’s surface elevation between 54° south and 60° north. The radar signal was reflected off the first object encountered, and was therefore affected by all constructed or natural objects on the Earth’s surface. For this reason, the SRTM data could be categorized as a digital surface model (DSM).
The SRTM Version 2 data were released in 2005, with certain artefacts fixed, including single-pixel errors and the noise that occurs in radar returns from water surfaces. About 3 arc-second SRTM data are available for most of the world’s continents, but 1 arc-second SRTM data are publicly available only for the United States. The specifications of SRTM data are provided in .
Table 1. The published profiles of the two data sets.
In this study, the NED data were used as the reference data to which SRTM would be compared. NED data are generated from the bare-Earth contours of USGS 1:24,000-scale topographic maps. The contours are then gridded at a constant spacing into quadrangle-based DEM tiles (Gesch et al. Citation2002; Gesch Citation2007). Each tile is clipped to the actual coverage extent (1° by 1°) to remove edge artefacts. The NED data are available at 1 arc-second, 1/3 arc-second and 1/9 arc-second postings on USGS’s seamless data warehouse.
Basing DEM points on equal intervals of latitude and longitude necessarily lead to east-west spacings that are only 83% of north-south spacings at the latitude of Los Angeles. For the purposes of this study, we reprojected both data sets to UTM Zone 11, resampling both to a spacing of 27.24 m to ensure that the DEM points are equidistant in both directions. In principle, the postings of the two data sets should coincide, though from the available documentation of these data sets it is clear that the implicit convolution functions are different: NED data are more point-like, with narrow convolution functions, while the SRTM elevations compare better to averages over pixel areas because of the wider implicit convolutions.
Before comparing the two data sets, it is important to point out that they use slightly different vertical data: NAVD88 for NED and EGM96 for SRTM. However, the differences between them are virtually constant over the study areas and as such would have no effect on the optimization of our approach; and they are submeter in magnitude. We therefore chose to ignore the difference.
5. Methods
5.1. Problem formulation
This section presents a formal mathematical representation of the components (slope, aspect, elevation, shift offset and shift angle) that influence the multi-source DEM differencing process. Our analysis is another form of that of Nuth and Kääb (Citation2011), with the exception of a constant term, which we didn’t consider since it has no influence on the optimization. Consider the same point 
in both SRTM and NED DEMs. The shift from SRTM to NED is a vector
, represented by its offset
and direction
, as Figure 3 shows. Let
denote the difference in elevations provided at a point by the two sources, and let
denote the error that would be observed at that point due to a misregistration of amount
alone. When
, the shift is in the direction of steepest slope, and because of misregistration the elevations obtained from the two sources should differ by
, plus any error due to mismeasurement. When
, the apparent slope parallel to
is
, therefore, the misregistration error caused by the shift equals:
(1)
Published online:
10 September 2015Figure 3. Vector representation of local shift 
in three-dimensional space.
-axis toward North,
-axis toward the East and z-axis measures the elevation above ground.
is the offset of the shift and
is the direction of the shift measured clockwise from North. The terrain slope is represented by an inclination
to the horizontal, and aspect is represented by an angle
measured clockwise from North.
Figure 3. Vector representation of local shift in three-dimensional space.
-axis toward North,
-axis toward the East and z-axis measures the elevation above ground.
is the offset of the shift and
is the direction of the shift measured clockwise from North. The terrain slope is represented by an inclination
to the horizontal, and aspect is represented by an angle
measured clockwise from North.
Our objective is to identify the unknowns d and β.
We assume that d and β vary slowly, and can be estimated using a block strategy. For example, a block of 1 km2 contains about 900 observations, and a block of 1.4 km2 contains about 2200 observations. Writing the problem as an optimization, and adding subscripts i to identify individual observations, the objective can be formalized as 
(2)
where n is the number of observations that fall inside the desired scene. Note that we assume a local shift for each block instead of a single global shift because the former emphasizes the effects of local geometric features (Flusser Citation1992; Goshtasby Citation1987), though we expect some degree of consistency in the estimates for each block. In the next section, the methods used for the shift estimation will be discussed.
5.2. Multidimensional unconstrained non-linear minimization
The following characteristics of the objective function indicate that this problem can be categorized as a multidimensional unconstrained non-linear optimization problem. First, the minimization of the objective function is decided by more than one variable (
and
), meaning that the problem space is multidimensional. Second, the objective function is non-linear since it is twice continuously differentiable. Third, the domain of permissible solutions is not limited to variables satisfying some equality or inequality constraints, therefore it is an unconstrained problem.
In the literature, there are many methods and algorithms developed to solve the multidimensional unconstrained non-linear problem, such as Newton’s method, Broyden’s method (Broyden Citation1965), line-search methods, trust-region methods (Celis, Dennis, and Tapia Citation1985), the Nelder–Mead simplex method (Nelder and Mead Citation1965), the conjugate-gradient method (Hestenes and Stiefel Citation1952) and their variants (Fan and Zahara Citation2007; Birgin and Martinez Citation2001). Most of these algorithms involve an iterative process that converges on a sufficiently accurate solution. Newton’s method is a typical example of the above processes. However, the method is very computational intensive [O(n3) operations per iteration] and requires many function evaluations at each iteration. Broyden’s method does not require the error-prone computation in Newton’s method and it is able to reduce the computational complexity from O(n3) to O(n2). But it may fail to converge on the optimal solution if poor starting points are selected. In comparison, line-search methods and trust-region methods introduce strategies, such as steepest descent (Wolfe Citation1969), to prevent the process from getting stuck at a local optimum, while preserving the advantages of local convergence in Newton’s and Broyden’s methods. All of these methods are derivative-dependent and are therefore computationally intensive. To overcome this drawback, the direct-search methods are built on sound heuristics and are derivative-free. Methods such as the Nelder–Mead algorithm and the conjugate-gradient method have drawn a lot of attention in the optimization field (Dennis and Schnabel Citation1989; Lewis, Torczon, and Trosset Citation2000). The conjugate-direction method is more suitable for solving problems with large numbers of variables, while the Nelder–Mead algorithm is very effective in solving minimization problems with few dimensions.
The Nelder–Mead algorithm begins by selecting a simplex of n + 1 points (n is the problem dimensionality, in our case, n = 2). Mathematically, the selection of n + 1 points is based on the fact that the derivative of the objective function of n variables could be estimated by the n + 1 function values through finite differences. There are four operations: (1) reflection, (2) expansion, (3) contraction and (4) shrink, for moving the vertices of the simplex until it converges. Reflection moves the location of the worst vertex through the centroid to its opposite side; worst means that the value of the objective function at that vertex is the least desirable. In our minimization problem, the worst vertex is the one with the largest objective value. By reflection, the solution space (the interior of the simplex) moves away from the worse set of values and moves closer to a better set. Expansion doubles the distance from the reflection point to the centroid. Similarly, contraction halves the distance from the reflection point to the centroid. Both the expansion and reflection operations reposition the reflected vertex and thus accelerate the search. If no more acceptable move of the vertex can be identified, the edge facing the best vertex is halved and the two vertices on the same edge move toward the best vertex to shrink and close up the simplex. Once the simplex converges to a stationary point, the optimal solution is found.
The unique features of Nelder–Mead allow it to discover patterns which cannot be immediately obtained from the original specifications, and the simplicity of this method can avoid the pitfalls of other sophisticated approaches. Therefore, it is selected as the algorithm in our study.
5.3. Workflow
Figure 4 illustrates the workflow for the DEM registration. This process is divided into three phases: data preparation, data processing, and results visualization:
Data preparation. This phase generates all required input data (slope, aspect and observed error) for the minimization process. First, the tiled SRTM and NED data were merged separately into a single data set. Then, the DEMs falling inside the two rectangular study areas, SGM and SMM, were extracted. As the NED data is assumed more accurate, the local slopes and aspects were estimated from NED data using the standard Horn (Citation1981) algorithm. Besides slope and aspect, the observed errors between SRTM and NED data sets were obtained by subtracting terrain height values on a pixel-by-pixel basis. Finally, these three images were converted to two-dimensional ASCII matrices as the input for the minimization process. ArcMap and its python library ArcPy (raster analysis module) are used to generate the terrain variables, including slope, aspect and the observed elevation errors.
Data processing. This phase optimizes the function introduced in Section 5.1 utilizing the Nelder–Mead algorithm in Matlab. As the premise of this procedure is that misregistrations are uniform in local blocks, the original data matrices were divided up into several blocks, and the shift computed separately in each block. Taking the 1 km × 1 km block as an example, and considering that we are using 30 m DEM data, then Block 1 (right matrix in Figure 5) will contain cells falling in the first 30 rows and 30 columns. Applying the shift model (in Section 5.1) at each pixel and summing over the 900 pixels, we obtained the function to solve for Block 1. Similarly, functions for other blocks can be obtained and the Nelder–Mead minimization can be applied to obtain one local shift vector per block.
Results visualization. Through the above process, we obtained the misregistrations between SRTM and NED for all blocks in terms of shift offset and angle. To visualize the shift in each block on the actual terrain, each such shift vector was added into a polyline shapefile. The start point of each shift vector is the centre of the scene. The end point of the shift vector is calculated from the offset and angle of the shift. As subpixel shift is usually small in comparison to the extent of the scene, the offset was lengthened by 25 times when being rendered on the map.
Published online:
10 September 2015
Figure 4. Three-phase process for DEM registration.
Published online:
10 September 2015Figure 5. Mapping from original data matrix into matrices of smaller blocks 1 km on a side.
Figure 5. Mapping from original data matrix into matrices of smaller blocks 1 km on a side.
6. Results
6.1. Experiment (A): method validation by synthetic data
To validate the proposed method, we constructed a cone-shaped mountain with steep terrain (the slope is 100%) to simulate similar steep terrain in our study area. The mathematical expression of the mountain is
subject to: 
where (0,0) is the centroid of the mountain base (a circle) and 
is the highest elevation of the synthetic mountain. The characteristic of the cone-shaped mountain is that the magnitude of slope is constant but the aspect is always changing. shows the spatial profile of the synthetic mountain. After constructing the mountain, a misregistration shift was introduced at 0.5 pixel in the eight directions of north, south, east, west, northeast, southeast, southwest and northwest (Figure 6 shows the shift from original location (coloured) to the northwest direction (grey)). Then the observed elevation difference, slope and aspect at each pixel cell were computed to feed into the Nelder–Mead algorithm for generating the estimated shift. Finally, the estimated shift and the actual shift were compared, giving the results shown in .
Table 2. Parameters of cone-shaped mountain.
Table 3. Comparison of estimated result by minimization and actual shift.
Published online:
10 September 2015
Figure 6. Manual shift of the synthetic mountain.
compares the result obtained from the Nelder–Mead estimation and the actual shifts for this mountain. The local shift vectors returned from the Nelder–Mead estimation were averaged to obtain the estimated shift. From , it can be seen that the proposed Nelder–Mead estimation shows good performance in identifying the shift in eight directions at a subpixel (½ pixel) level. Errors in the estimated offset were no more than 2% and errors in shift angle less than 1°. This experiment proves the validity of the proposed method in handling the co-registration of DEMs with terrains on which aspect changes frequently. We also plotted (Figure 7) the elevation differences averaged over two aspects (east and northwest) when the data has misregistration in cardinal (e.g. east) and ordinal directions (e.g. northwest). These sinusoidal curves indicate that the averaged elevation differences maximize at the aspect (direction) of introduced shift. This conclusion can be deduced from the spatial model introduced in Section 5.1 – at the direction of the shift, the elevation difference is toward the steepest slope and therefore, the elevation difference at that direction is maximal. This finding can be used to assess the accuracy of estimated shift from the real data sets.
Published online:
10 September 2015Figure 7. Averages of elevation difference caused by horizontal shift in terms of aspect: (a) represents case of shift toward east and (b) represents case of shift toward northwest.
Figure 7. Averages of elevation difference caused by horizontal shift in terms of aspect: (a) represents case of shift toward east and (b) represents case of shift toward northwest.
6.2. Experiment (B): local shift identification in both SGM and SMM regions
Experiment (A) testified to the validity of the proposed method in identifying subpixel misregistration. This experiment explores the actual shifts that exist between SRTM and NED data in both study sites: SGM (Figure 8(a)) and SMM (Figure 8(b)). The block size for both regions is 1.4 km × 1.4 km. This block size was chosen to illustrate the misregistration because using this block size, the result obtained is the most optimal (has the highest objective values), in comparison to other block sizes such as 1 km. The numbers on the arrows refer to the offsets of the shifts (in meters) and the directions refer to the angle of the shifts. From the figures, we can tell that the SRTM data set maintains a consistent subpixel shift toward the northwest compared to the NED data set within both regions. We also found that the offset (25.8 m) of the shift in the SGM region, although also subpixel, is greater than that in SMM region, which has less steep terrain.
Published online:
10 September 2015Figure 8. Visualization of identified shift at 1.4 km × 1.4 km scale from both SGM (a) and SMM (b) regions.
Figure 8. Visualization of identified shift at 1.4 km × 1.4 km scale from both SGM (a) and SMM (b) regions.
On average, the shifts are (25.8 m, −23.6°) for the SGM region and (17.89 m, −27.6°) for the SMM region at the tested scale. By using blocks rather than estimating a single shift for each region, we are able to explore the possibility of spatial variations in shifts. We can infer from Figure 8 and the above experiments that the misregistration is a systematic offset in this case, in two continuous terrains that have similar topographical features. However, the misregistration is not a global offset; its value tends to be larger in steeper terrains (such as the SGM region) and smaller in less steep terrains (such as the SMM region). Figure 8(a) and 8(b) shows possibly interesting patterns of variation between the shifts estimated from blocks.
One might argue that the selection of block size in the experiment may affect the results. In order to explore a possible effect, the optimization procedure was performed on 14 different block sizes with linear dimensions: 250 m, 500 m, 1 km, 1.2 km, 1.4 km, 1.6 km. 1.8 km, 2 km, 3 km, 4 km, 5 km, 6 km, and 7 km; and study area as a whole. We found that the shifts obtained from different block sizes are all similar, with offset falling in ranges between 24.3 m and 26.6 m and angle ranges from −24.5° to −23.2° in the SGM region. The misregistrations obtained in the SMM region fall between 16.78 m and 19.13 m of offset and the angle of the shift ranges from −28.78° to −26.02°.
Besides the multi-scale analysis, we also generated a map of the elevation differences between the original SRTM and NED data. Figure 9 shows these differences averaged and plotted by aspect for 1.4 km blocks. The peak differences for the SGM region occur when aspect equals 337° and for the SMM region when aspect equals 332° These values are consistent with the shifts (−23.6° in the SGM region and −27.6° in the SMM region) identified during the registration procedure. This experiment further verifies the effectiveness of the proposed method in detecting the misregistration.
Published online:
10 September 2015Figure 9. Average of elevation differences between SRTM and NED plotted against aspect in the SGM (a) and SMM (b) regions.
Figure 9. Average of elevation differences between SRTM and NED plotted against aspect in the SGM (a) and SMM (b) regions.
7. Discussion and conclusion
This article presents a robust optimization technique to identify the spatial variation in horizontal misregistration between multi-source DEMs, through an extension of the method proposed by Nuth and Kääb (Citation2011). We introduced the formulation of the problem by proposing a block-specific shift vector that considers elevation differences, local slope and aspect. The Nelder–Mead algorithm applied to solve the multidimensional unconstrained non-linear minimization problem produced satisfying results through a number of experiments using both synthetic data, and real data sets in the SGM and SMM regions.
Quantification of horizontal misregistration on a subpixel level is essential for performing precise scientific analyses. From the multi-scale analysis, we estimated an average horizontal shift in the SGM region of 25.8 m with 23.6° deviation to the north; for the SMM region, the average horizontal shift is 17.89 m with 27.6° deviation to the north. We also found that the horizontal misregistration is a relatively systematic offset in continuous terrains that have similar topographical features. However, the misregistration is not a global offset; its value tends to be larger in steeper terrains (such as the SGM region) and smaller in less steep terrains (such as the SMM region).
The major contribution of this article is thus the extension of the method of Nuth and Kääb (Citation2011) to blocks, allowing the spatial variation in offset to be estimated. We found that such variation can be estimated satisfactorily in hilly terrains. Where the terrain is relatively flat, the error term estimated in the model tends to be close to 0 (affected by the 
term in Equation 1), making it hard to detect the misregistration pattern during the optimization simulation. Therefore, it would be difficult to estimate offsets using this method. This raises an issue for implementation, when corrections are applied to real data. If sharp changes in estimated offsets occur at block boundaries, there is the danger that this will introduce unwanted effects in the resultant DEMs. Thus for implementation, we suggest that estimates for relatively flat blocks be replaced by interpolated estimates from nearby hilly blocks. Clearly this will be problematic in large areas of flat terrain, but in such areas DEM misregistration is less likely to cause problems in applications.
Technically, one difficulty to overcome in solving the optimization problem is local optimality due to the non-convex nature of this problem. During the experiment, we found some evidence of local optima in some blocks. For example, the original shift identified in the SGM region for the block at row 3 and column 11 was (34.47 m, −23.23°), which implies an overestimation of the offset, given that the averaged shift for the SGM region is (25.8 m, −23.6°). When such indications of a local optimum occurred, the optimization algorithm was rerun on the block by using the average shifts as an initial state. If the new estimates led to a smaller objective function, that is, a better solution, the original local estimates would be replaced. In this way, at least some of the local optima could be avoided.
However, after removing the local minima, there remain some anomalies in the shifts. For instance, consider the obviously abnormal shifts identified in two scenes numbered 1 and 2 in the black box in Figure 8(b). By looking into the actual satellite images of both blocks, we found that substantial portions of both blocks are developed, and there are large warehouses in block 1 and many residential buildings in block 2. The elevation surface containing the heights of the buildings in the SRTM data is not consistent with the actual terrain surface anymore. Development and other inherent uncertainties in DEM data sets, such as those caused by the coarse resolution, are probably the reason for these anomalies, along with the general flatness of the terrain in these blocks. Nevertheless, using the proposed method, we are still able to produce valid and satisfying results for most of the regions in our study area.
In the future, we will work on quantifying uncertainties (such as terrain complexity and existence of build-ups and non-natural features) that cause the anomalies in the estimation of DEM misregistration, in order to produce more accurate results. Then we will expand this research to detect misregistrations in urban areas. Both of the research directions require an effective mechanism in detecting the locations of human construction, and a method to remove the effects caused by these features in the DSM data. Land-use data indicating the footprints of buildings need to be considered and a more complicated optimization procedure needs to be developed. We will also work on improving the computing performance of the optimization process. For the study areas used in this article, a single workstation (Intel Core i5 CPU with 2*3.2 GHz cores; 8 G memory) can finish processing the DEM registration using both SRTM and NED data set in about 0.2 h. When dealing with the problem over a much larger extent and with finer-resolution data, a single workstation will be far less satisfying in computing performance. As the optimization is conducted on equally sized regions, it fits naturally the divide-and-conquer paradigm employed in GPGPU (general purpose graphics processing unit) processing. Therefore, parallelization and optimization of the algorithm, and its adaptation to GPU clusters, will be another focus of our further work.
Disclosure statement
No potential conflict of interest was reported by the authors.
| DEM source | NED | SRTM |
|---|---|---|
| Full name | National Elevation Data | Shuttle Radar Topography Mission |
| Provider | USGS | NASA |
| Type | Bare-earth elevation | Surface elevation |
| Spatial resolution | 1/9 arc-second (~3 m) 1/3 arc-second (~10 m) 1 arc-second (~30 m) | 1 arc-second (~30 m) 3 arc-second (~90 m) |
| Vertical accuracy | 2.44 m (RMSE) | 9.0 m (90% confidence level) (Rodríguez, Morris, and Belz Citation2006) |
| Horizontal accuracy | n/a | ±20 m (90% confidence level) (SRTM Mission Statistics Citation2000) |
| Parameter | Notation | Value |
|---|---|---|
| Height |  | 2 |
| Slope |  | 100% |
| Aspect |  |  |
| Pixel size |  | 0.01 |
| Relative shift |  | (2.732 pixel size, 30°) |
| Absolute shift |  | (0.02732, 30°) |
| Number of blocks | 1 × 1 |
| Case | Shift direction | Shift distance | Actual shift | Estimated shift |
|---|---|---|---|---|
| (1) | North | ½ pixel | (0.005, 0°) | (0.005, 0.1°) |
| (2) | South | ½ pixel | (0.005, 180°) | (0.005, 180.2°) |
| (3) | East | ½ pixel | (0.005, 90°) | (0.005, 90°) |
| (4) | West | ½ pixel | (0.005, 270°) | (0.005, 269.83°) |
| (5) | Northeast | ½ pixel | (0.005, 45°) | (0.0049, 45.78°) |
| (6) | Southeast | ½ pixel | (0.005, 135°) | (0.005, 135.02°) |
| (7) | Southwest | ½ pixel | (0.005, 225°) | (0.0049, 225.97 °) |
| (8) | Northwest | ½ pixel | (0.005, 315°) | (0.005, 315.08°) |