ESA Earth Home Missions Data Products Resources Applications
EO Data Access
How to Apply
How to Access
FRINGE '96 Workshop: ERS SAR Interferometry, 30 September - 2 October 1996
Site Map
Frequently asked questions
Terms of use
Contact us



Multiresolution Analysis of DEMs: Error and Artifact Characterization

Dragos Luca Image Science Division, Communication Technology Laboratory, ETHZ, Gloriastr. 35, CH-8092 Zürich, Switzerland
Mihai Datcu Deutsche Forschungsanstalt für Luft und Raumfahrt e.V., Deutsches Fernerkundungsdatenzentrum, D-82234 Wessling, Germany
Klaus Seidel Image Science Division, Communication Technology Laboratory, ETHZ, Gloriastr. 35, CH-8092 Zuerich, Switzerland


Digital elevation data is often affected by errors due to excessive interpolation or to interferometric artifacts. The local analysis of the DEM roughness allows both the detection of these errors and the segmentation of the elevation data for a better understanding of the existing geological structures. This analysis can be performed by means of fractal dimension estimators. We compare several fractal estimation methods and show that those based on the multiresolution (wavelet) data analysis yield the best results from the point of view of their segmentation capabilities. This is a natural conclusion considering that fractals and wavelets rely on many common concepts. As an example we show an application in which SAR intereferometric artifacts are revealed and the elevation data is separated in different roughness classes using fractal dimension measurements and an unsupervised clustering algorithm.
Keywords: DEM, fractals, wavelets, multiresolution, artifacts, segmentation

1 Introduction

Digital Elevation Models (DEMs) have become an important tool in many remote sensing applications like e.g. SAR simulators, orthorectification of satellite and airborne images, classification of ground cover types, Geographic Information Systems, etc. Still in most of the cases, the resolution of the DEMs derived from digitized topographic maps is unsufficient for the requirements of the corresponding applications. SAR interferometry allows the achievement of much better resolutions, however this technique is not yet available on a large scale and the resulting surface is affected by typical artifacts. The detection of the artifacts is an important task, especially in remote areas where no other digital elevation data is available for comparison.

Relying on the observation that relief conserves the same statistical characteristics over a wide range of scales, some DEM characterization techniques have been developed that perform a multiresolution analysis of the elevation data. These techniques use fractal models to measure the local roughness of the DEM and extract derived information. In this way, better interpolation methods for higher resolution elevation data sets have been developed (Yokoya and Yamamoto 1989, Polidori et. al. 1991, Franceschetti et. al. 1994) and a better understanding of the geological and geomorphic processes has been achieved (Huang and Turcotte 1990,Clarke and Schweizer 1991,Clarke 1988).

This work investigates the use of multiresolution techniques for the characterization of DEM roughness and points out a method for the detection of errors and artifacts in intereferometric elevation models. The roughness of a DEM is estimated locally by measuring the fractal dimension of the underlying model. Two fractal dimension estimators are compared: the power spectrum estimator and the wavelet-based estimator. The power spectrum method was selected as a reference since it has often been designated to be the algorithm which achieves the best performance (Stewart et. al. 1993, Schepers et. al. 1992). Our approach concentrates on the wavelet-based method and is motivated by the fact that wavelets and fractals are closely related by the concept of scale and share many common properties (Akujuobi and Baraniecki 1994). We compare the two methods both on synthetic and on real elevation data. The result of this comparison shows that the wavelet-based estimator achieves a better reliability of the measurements in terms of their standard deviation and is better suited for the segmentation of fractal images.

The paper is organized as follows. Section 2 gives an overwiev of the fractal and wavelet theory emphasizing the common concept of scale. Section 3 presents the results of a comparison between the spectral and the wavelet-based fractal estimators based on synthetic images. This comparison is performed in terms of image segmentation capabilities, i.e. the algorithms are applied to image windows of small size and the obtained set of measurements is analysed statistically. In the last section we discuss some examples in which the fractal analysis of DEMs reveals artifacts in the computation of elevation data and allows the separation of different roughness classes.

2 Elements of Fractal and Wavelet Theory

Fractals are mathematical objects that show the same structure when examined at all possible scales (Mandelbrot 1982). Although more formal definitions can be given, this qualitative characterization expresses the very essence of the fractal phenomenon and at the same time represents the basis for all fractal analysis algorithms: the algorithms check for fractal behaviour simply by examining the object at several scales (resolutions). The basic parameter characterizing a fractal object is its dimension: while non-fractal objects have dimensions given by integers (1 for curves, 2 for surfaces, etc.), fractals will have fractional dimensions since they represent transition structures between curves and surfaces, surfaces and solid bodies, etc.

The statistics of roughness measurements has been shown to agree - in a limited resolution range - with that of specific fractal models and several attempts have been made to characterize DEMs by means of the fractal theory. The most popular model used in this respect is the fractional Brownian motion (fBm) model (Voss 1988). This model describes a signal B(t) characterized by the fact that its increments between two moments of time t1 and t2 have a variance proportional to a power 2H of the time lag |t2-t1|. The parameter H is called "Hurst exponent" (0<H<1) and is related to the dimension D of the corresponding fractal by (Voss 1988) D = n + 1 - H, n being the topological dimension of the space in which the fractal object is represented (n=1 for fBm time-functions, n=2 for fBm surfaces, etc.). D will typically measure the "roughness" of the signal and most fractal analysis algorithms concentrate on an accurate estimation of D. Note that fBm signals are nonstationary; this makes their analysis both from the theoretic and from the practical point of view quite difficult. Several methods have been developed that try to cope with these difficulties.

A classical algorithm for the estimation of D is based on the computation of the time (space) averaged power spectrum of the signal. This power spectrum is represented in a log-log plot vs. frequency. It can be shown (Voss 1988, Saupe 1988) that if the signal is a fBm the points of the plot will align along a straight line. The slope of this line is directly related to the dimension D of the fractal and is usually estimated via a regression technique. This method of fractal dimension estimation is called the power spectrum method and has been shown to yield very accurate estimation values (Stewart et. al. 1993).

But other methods exist as well and deserve increasing attention since they present some very useful features for different applications. Our research focuses on an alternative way for the estimation of D based on the wavelet decomposition of fBm signals. Just like fractals, wavelets heavily rely on the scale concept, i.e. on the analysis of the data at several resolutions (Mallat 1989, Wornell 1993). Given a signal x(t), a sequence A[m]x(t) of approximations to x(t) is constructed, each A[m]x(t) representing the signal approximation at a given resolution (scale) m. If we define D[m]x(t) the detail signal at resolution m to be the difference between the two succesive approximations A[m+1]x(t) and A[m]x(t), then this detail signal can be written as an orthogonal expansion of the signal x(t) using some special basis functions that are called "wavelets". The wavelet transform represents thus a way of quantifying signal changes from one scale to another. The wavelets themselves are obtained using the concept of scale: they are all dilations and translations of a single function called "mother wavelet".

The similar fundamental concepts of the fractal and wavelet theory have stimulated several researchers to investigate this relation in more detail. The main result in this field is due to G. W. Wornell who showed that the wavelet transform applied to a fBm signal whitens the signal, i.e. the transformed-domain samples at a given resolution m (the detail signal samples) become stationary and weakly correlated and both the theoretic and the numeric analysis are much easier to perform. The direct estimation of the fractal dimension relies in this case on the the computation of the detail signal variances at different scales. A maximum likelihood estimation algorithm (Wornell and Oppenheim 1992) can be used to estimate from several variance measurements the fractal dimension D.

3 Comparisons of Fractal Models

In order to compare the performance of the methods presented in the previous section, a set of 256-by-256 spatially isotropic fBm surfaces ranging from dimension 2.0 to dimension 3.0 was generated. In each image, the dimension was estimated locally in a sliding window of size 32-by-32 pixels with the spectral and the wavelet-based algorithm (using Daubechies wavelets with 4 filter coefficients).

The comparison of the different estimation algorithms is performed mainly by means of two parameters: the "mean of measurements" and the "uniformity of measurements". These parameters are plotted vs. the real dimension of the fractal surface in figure 1. The mean of the measurements represents the mean value of the fractal dimension estimations for all the positions of the sliding window. For a uniform surface of fractal dimension D, this mean should yield exactly the value D, i.e. the ideal shape of the mean plot is that of the diagonal line. The uniformity of measurements is defined as the standard deviation of the measurements about the measured mean. The uniformity plot should thus be flat and of minimal value to characterize a reliable estimation method.

The results of this analysis show that although the spectral method has a better performance in terms of measurement accuracy, the wavelet-based method has a substantially better uniformity. A good uniformity is in our opinion more important for image segmentation than the exact estimation of the numeric values for the fractal dimension, provided that the relation mean vs. D remains monotonic. This is the case for the wavelet-based method, where the mean vs. D plot can be used as a LUT for the correction of estimated values.

Figure 1: Mean and uniformity of fractal dimension measurements for synthetic surfaces. Window size is 32 x 32 = 1024 pixels

These results can be visualized by means of an example. Figure 2a shows a synthetic DEM (shown as a shaded relief) consisting of a central square of fractal dimension D=2.8 surrounded by a rectangular border of fractal dimension D=2.2. Figures 2b and c show the estimation results for D (i.e. for the roughness of the relief) using the spectral and the wavelet-based algorithms respectively in a sliding window of size 32 x 32 pixels. The estimated numeric values are scaled linearly from [2.0, 3.0] to [0, 255] and are presented as gray scale images. White corresponds to a high fractal dimension ("rough") and black to a low fractal dimension ("smooth"). To avoid problems generated by discontinuities on the borders of the image, we did not compute the fractal dimension in a swath of 16 pixels on each image side. This swath is shown here in black.

Figure 2: Segmentation example for synthetic DEM consisting of a central square of fractal dimension D=2.8 surrounded by a border of fractal dimension D=2.2. The dimension is estimated in a sliding window of size 32 x 32 pixels. a) Original image, b) Fractal dimension estimated with the spectral algorithm, c) Fractal dimension estimated with the wavelet-based algorithm

Figure 2 confirms the conclusions of the mean and uniformity plots. Due to the better uniformity of the wavelet-based estimation algorithm the fractal dimension map in figure 2b appears less noisy than the one computed with the spectral method (figure 2c). This leads to a better visual appearance and to an easier thresholding for segmentation.

4 Application to Digital Elevation Models

The algorithms described in the previous section have been tested on two elevation data sets that present different characteristics and allow us to consider different aspects of the fractal estimation process.

The first data set consists of three DEMs of the region of Davos, Switzerland at resolutions of 50m, 25m and 10m respectively (figure 3a). The DEMs were obtained by digitization of elevation data from maps of different resolution. For each of the DEMs the fractal dimension was computed in a sliding window of a size selected roughly proportional to the dimensions of the DEM (16 x 16 pixels for the 50m DEM, 32 x 32 pixels for the 25m DEM and 64 x 64 pixels for the 10m DEM).

Figure 3b shows the histograms of the estimated values for D obtained using the wavelet-based method. First note that all the histograms show a single peak. This fact means either that the whole considered area has a single fractal dimension or that it consists of several areas with very close fractal dimensions which cannot be differentiated by this method. Note also that the peak is at about the same position for all three histograms, i.e. the structure of the terrain is invariant to scale. Since the data of the three DEMs was obtained by independent processes and not by interpolation of one set to another, this example proves that the terrain shows indeed a fractal structure. The dimension of this fractal lies in the range of D=2.1 to D=2.2. Similar fractal dimension values have been reported for DEMs of other areas as well (Polidori et. al. 1991, Clarke and Schweizer 1991).

Figure 3c shows the histogram of the estimated values for D using the spectral method. The variance of the measurements is higher and the estimation of D is less reliable. This result confirms the simulations described in the previous section.

Figure 3: Fractal dimension measurements for three DEMs of different spatial resolution of the Davos area. a) DEMs presented as a shaded relief: left to right 50m, 25m, 10m; b) Histograms of fractal dimension measurements for the DEMs in figure a) using the wavelet-based method? c) Histograms of fractal dimension measurements for the DEMs in figure a) using the spectral method

Our second example shows the possibility to use the estimation of fractal dimensions to detect artifacts in elevation data. The data consists of two DEMs of an area along the river Rhine in Germany. The first DEM (figure 4a) was obtained from digitized cartographic information and has a resolution of 25m, while the second one (figure 4b) was obtained by SAR interferometry and has a resolution of about 20m. Note also that the topographic DEM is represented in geographical coordinates and the interferometric DEM is shown in range/azimuth coordinates which are slightly rotated with respect to the geographical ones. In spite of these differences the data was not preprocessed (e.g. reinterpolated and rotated) since that would have affected the fractal structure of the images.

Figure 4: a) Topographic DEM of an area along the river Rhine; b) Interferometric DEM of approximately the same area

The topographic and the interferometric DEM are analysed by computing the fractal dimension in a sliding window of size 32 x 32 pixels. The histograms of the values for D obtained with the wavelet method are shown in figure 5a and b respectively. While for the topographic DEM, D lies in the range 2.0 ... 2.4 (the "usual" values for the fractal dimension of terrain structures), the interferometric data has much higher fractal dimensions, ranging from 2.4 to 2.8. These values indicate that the interferometric method induces some artifacts leading to an unusually high roughness of the elevation data. The estimation of the fractal dimensions provides thus a good method to detect the artifacts in an automatic way.

Figure 5: a) Histogram of fractal dimension measurements for the topographic DEM; b) Histogram of fractal dimension measurements for the interferometric DEM

Additionally, the histograms of D allow the segmentation of the DEMs according to their roughness. An unsupervised classification algorithm developed by Narendra and Goldberg (Narendra and Goldberg 1977) looks first for the dominant modi of the image histogram, then clusters the pixels around these modi, setting thresholds in the valleys between them. The algorithm was applied for the histograms in figure 5 and resulted in the classification map of figure 6. As expected, for both the topographic and the interferometric DEM, the clustering algorithms detected 2 classes of roughness: the first one corresponds to the mountain area, the second one to the plain area. These classes are presented here in pseudocolors to allow for a better visual separation.

Figure 6: a) Unsupervised classification of the topographic DEM according to its roughness, b) Unsupervised classification of the interferometric DEM according to its roughness

As a final experiment, another roughness measure, the standard deviation of the elevation values in an image window is computed for the topographic and the interferometric DEMs and compared to the fractal dimension. The result of this comparison is shown as a scatterogram of the values of the standard deviation vs. D (figure 7). Little correlation can be seen between the two roughness measures, especially in the case of the interferometric DEM where the points are scattered all over the plane. Obviously, the fractal dimension characterizes the roughness of the relief in a different way than the standard deviation of the elevation values. While the standard variation is a parameter that assumes a Gaussian distribution, D relies on a fractional Brownian motion model which is non-Gaussian. The appropriateness of one of these measures is related to the assumptions that can be made on the statistics of the image. Usually these statistics are not checked and an a priori assumption is made, based on some physical and/or empiric considerations. We feel that in this context a more rigorous approach is necessary, implying some better mathematical modeling. This approach is the object of further study and research to be reported soon.

Figure 7: Scatterogram of the standard deviation of elevation data vs. fractal dimension for the topographic DEM b) Scatterogram of the standard deviation of elevation data vs. fractal dimension for the interferometric DEM

5 Conclusions

The present paper presents an automatic method for the detection of errors and artifacts in the generation of interferometric DEMs. The detection is based on the estimation of terrain roughness and uses the observation that zones where interferometric artifacts appear have a higher roughness than "normal" terrain.

The roughness is measured assuming a fractal model for the elevation data and is expressed by the fractal dimension, estimated locally in a small sized window. Several estimation methods are compared from the point of view of their segmentation capabilities. Based on important conceptual similarities between fractals and wavelets, we make the hypothesis that wavelets are a very appropriate analysis tools for fractal images. This assumption is confirmed by experiments both on simulated and on real elevation data: wavelet-based fractal dimension estimators show the highest reliability in terms of measurement variances and are thus better suited for DEM characterization than other methods.


We would like to acknowledge the German Remote Sensing Data Center (DFD) at DLR, Wessling, Germany for providing us the DEMs of the Rhine valley and the Remote Sensing Laboratories of the University of Zürich, Switzerland for providing us the 10m resolution DEM of the Davos area. Special thanks also to Dr. Marcus Schwäbisch from DLR for the derivation of the interferometric SAR DEM.


Akujuobi, C. M. and Baraniecki, A. Z. (1994).
A Comparative analysis of wavelets and fractals. In Leondes, C. T., editor, Digital Image Processing: Techniques and Applications, volume 67 of Control and dynamic systems, pages 143-197. Academic Press, Inc.
Clarke, K. C. (1988).
Scale-based simulation of topographic relief. The American Cartographer, 15(2):173-181.
Clarke, K. C. and Schweizer, D. M. (1991).
Measuring the fractal dimension of natural surfaces using a robust fractal estimator. Cartography and Geographic Information Systems, 18(1):37-47.
Franceschetti, G., Mihliaccio, M. and Riccio, D. (1994).
SAR simulation of natural landscapes. In stein, T. I., editor, Surface and Atmospheric Remote Sensing: Technologies, Data Analysis and Interpretation, IGARSS'94, Pasadena, pages 1181-1183.
Huang, J. and Turcotte, D. L. (1990).
Fractal image analysis: application to the topography of Oregon and synthetic images. J. Opt. Soc. Am. A, 7(6):1124-1130.
Mallat, S. G. (1989).
A theory for multiresolution signal decomposition: The wavelet representation. IEEE Trans. on PAMI, 11(7):674-693.
Mandelbrot, B. B. (1982).
The Fractal Geometry of Nature. Freeman, CA.
Narendra, P. M. and Goldberg, M. (1977).
A non-parametric clustering scheme for Landsat. Pattern Recognition, 9:207-215.
Polidori, L., Chorowicz, J. and Guillande, R. (1991).
Description of terrain as a fractal surface and application to digital elevation model quality assessment. Photogrammetric Engineering and Remote Sensing, 57(10):1329-1332.
Saupe, D. (1988).
Algorithms for random fractals. In H.-O. Peitgen and D. Saupe, editors, The Science of Fractal Images, chapter 2, pages 71-120. Springer-Verlag.
Schepers, H. E., van Beek, J. H. G. M. and Bassingtwaighte, J. B. (1992).
Four methods to estimate the fractal dimension from self-affine signals. IEEE Engineering in Medicine and Biology, pages 57-71.
Stewart, C. V., Moghaddam, B., Hintz, K. J. and Novak, L. M. (1993).
Fractional Brownian motion models for synthetic aperture radar imagery scene segmentation. Proceedings of the IEEE, 81(10):1511-1521.
Voss, R. F. (1988).
Fractals in nature: From characterization to simulation. In H.-O. Peitgen and D. Saupe, editors, The Science of Fractal Images, chapter 1, pages 21-70. Springer-Verlag.
Wornell, G. W. (1993).
Wavelet-based representations for the 1/f family of fractal processes. Proceedings of the IEEE, 81(10):1428-1450.
Wornell, G. W. and Oppenheim, A. V. (1992).
Estimation of fractal signals from noisy measurements using wavelets. IEEE Trans. on Signal Processing, 40(3):611-623.
Yokoya, N. and Yamamoto, K. (1989).
Fractal-based analysis and interpolation of 3D natural surface shapes and their application to terrain modeling. Computer Vision, Graphics and Image Processing, 46:284-302.

Keywords: ESA European Space Agency - Agence spatiale europeenne, observation de la terre, earth observation, satellite remote sensing, teledetection, geophysique, altimetrie, radar, chimique atmospherique, geophysics, altimetry, radar, atmospheric chemistry