Multiresolution Analysis of DEMs: Error and Artifact Characterization
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
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.
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
|Copyright 2000 - European Space Agency. All rights reserved.|