Introduction

1As a Digital Elevation Model (DEM) is a representation of a surface, several algorithms have been developed to study the surface properties of DEMs and to provide a large set of descriptor attributes (Wilson and Gallant, 2000). Given that altitudes correspond to grey levels, it is then possible to describe, quantify, and model rough surfaces using image-processing techniques in the spatial and frequency domains and to apply pattern recognition processes. Textural analysis is closely linked to roughness assessment. The concept of texture is quite difficult to define as it often includes the notions of roughness, regularity, contrast, and thinness. Numerous authors have attempted to outline and clarify this notion. J. C. Russ (1999) proposed that roughness equates with the high frequencies of the 2D signal constituting a DEM. He pointed out that a DEM contains three levels of information. The first level, related to the lower frequencies, coincides with its form; the second level, related to the middle frequencies, corresponds to its waviness; and the third level corresponds to its roughness.

2Some recent techniques in image analysis use fractal and/or multi-fractal approaches to characterize the texture of a greyscale image or the roughness of a surface. Concerning DEM surfaces, several authors have demonstrated that the topography of the Earth generally exhibits fractal characteristics and that the relief preserves the same statistical characteristics over a wide range of scales (Huang and Turcotte, 1989; Klinkenberg and Goodchild, 1992). By describing the terrain as a fractal surface, the local or the global analysis of DEM roughness reveals interpolation artefacts, provides an assessment of the DEM quality (Polidori et al., 1991; Datcu et al., 1996), and can assist in understanding erosional phenomena (Chase, 1992; Cheng et al., 1999).

3As every landscape appears to have a particular fractal dimension value, it implies that this dimension calculated locally has to be related with the particular features of the local landscape units. On the other hand, each local landscape unit presents its own morphologic feature in relation to the erosion processes and mainly to the nature of the basement. Thus, the purpose of this research is to study the relationships between geomorphic features and the surface roughness of a DEM by locally measuring fractal dimension in 3D space.

4The purpose here does not include demonstrating the fractal behaviour of the DEM nor comparing the performances of different fractal surface estimators. The first section of the paper presents a brief overview of fractal geometry and explains the notion of local fractal dimension. The second section describes the procedure, while the third section reports the results obtained in two case studies. These training areas have been chosen according to diverse lithological conditions. The well known geological features of the sedimentary Vittel region leads to investigate the ability of the treatment to extract these features in accordance to their morphology. Taking into account the former results, the treatments applied to the Mount Ararat illustrate the reliability of the method to emphasize local roughness differences related to regional tectonism and to detect peculiar and unknown features that can be useful for geological and geomorphological mapping.

Fractal geometry and local fractal dimension

Fractal geometry

5Fractal geometry, introduced and developed by B. Mandelbrot (1982), provides a mathematical description of a wide range of natural forms and phenomena. Fractal objects are defined as scale-invariant (self-similar or self-affine). This means that the fractal object can be presented as an assemblage of rescaled copies of itself. Self-similarity occurs when the rescaling is isotropic or uniform in all directions, and self-affinity occurs when the rescaling is either anisotropic or dependent on direction.

6Fractal objects exhibit details at arbitrarily small scales, and they are too complex to be represented in a Euclidean space. Also known as the Hausdorff-Besicovitch dimension (Falconer, 1990), the fractal dimension differs from the more familiar Cartesian or topological dimension. In this last case, integer values are required: 1 for a line, 2 for a surface, and 3 for a volume. The fractal dimension measures the complexity of the object. A shape with a higher fractal dimension is more complicated or irregular than one with a lower dimension. For example, a shape with a fractal dimension falling within the range [1, 2] fills more space than a one-dimensional curve and less space than a two-dimensional surface. Self-similarity is defined statistically when it cannot be tested through an infinite range of scales. The statistical fractal behaviour is then related to a given scale range. When the fractal dimension of a particular pattern changes within consecutive ranges of scale, one generally refers to the notion of multi-fractality.

7Several methods are available for estimating the fractal dimension of surfaces, such as the fractional Brownian model (Mark and Aronson, 1984), triangular prism areas (Clarke, 1986), box-counting (Falconer, 1990), and the projective covering method (Xie and Wang, 1999). Measurement may either be direct, when it is applied to grey tones, or indirect as in the case of the examination of “profiles” (one-dimensional transects) or isolines (conversion of the surface into a contour map). In this study, the box-counting method was used because it can be applied to various sets of any dimension and patterns with or without self-similarity (Peitgen et al., 1992). According to K. Falconer (1990) who has discussed the mathematical aspect of the box-counting method, the fractal dimension D can be derived from the relation:

81=N_ss^D or D—logN_s/log(s) where N_s corresponds to the number of boxes of size (s) needed to cover the structure. The D value is calculated with the following formula:

Local fractal dimension

9The local measurement of fractal dimensions has been developed to investigate image texture and surface roughness. When this type of measurement is applied to an image I, the results are reported in an image R where the value of each point (i,j) corresponds to its fractal dimension FD calculated in a window of size (2w + 1) x (2w + 1):

10Textural features derived from the local measurement of the fractal dimension have been used for segmentation and classification purposes. Among these techniques, the Fourier power spectrum (Pentland, 1984), the blanket method (Dellepiane et al., 1991), differential box-counting (Sarkar and Chaudhuri, 1994; Chaudhuri and Sarkar, 1995), and the fractional Brownian motion (Chen et al., 1989; Toennies and Schnabel, 1994) are the most common. Using a wavelet technique, M. Datcu et al. (1996) have estimated the local roughness of various DEMs and provided some examples showing that the fractal analysis of these DEMs allows different roughness classes to be distinguished and some artefacts, due to the computation of elevation data, to be detected. For computing the local self-similar properties of a DEM, Y.C. Cheng et al. (1999) have developed a 3D box-counting method applying the triangular prism surface method proposed by K.C. Clarke (1986). They observe that fractal dimension values vary as a function of altitude and have interpreted this phenomenon as reflecting spatial variability in erosional potential. This last method differs from the previous ones because the fractal dimension is estimated in unit areas and not at each point of the image. In our approach, the fractal dimension of each pixel of the DEM surface is calculated inside a moving window centred on this pixel, by using a 3D box-counting adaptive method. The advantage of this treatment is that the volume corresponding to DEM section observed locally is directly related to the altitude values of the pixels defining the DEM surface. The treatment does not require any modelling or interpolation of the surface in order to calculate this dimension.

Procedure

11Inside a cube of size s x s x s centred on the pixel of the studied DEM, the volume corresponding to the surface is defined by a set of voxels. A voxel is an elementary cube, the sides of which are equal to the pixel size. Thus, in the (x,y) space, each point (i,j) of the cube contains a number vs(i,j) of voxels falling between 0 and s. The total number of voxels describing the volume is equal to:

12It is easier to calculate vs (i,j) in a bi-dimensional space as follows:

13where I is the original image, Ic the value of the central pixel, ps the pixel size, h a coefficient that defines the vertical resolution, and s the cube size. The cube is partitioned into boxes of size q varying between 1 and s/2 (fig. 1). Each of these boxes is considered as filled if at least one voxel is contained in this box. In other words, the maximum of vs(i,j)/q in each cell is determined. The computation is done as follows:

14where Max is a function calculating the maximum. The fractal dimension corresponds to the inverse of the slope P=ln(q)/ln(Ns), where q is the size of the box and Ns the total number of filled boxes. When calculating the slope, the coefficient of determination, R², is computed. The estimate of the fractal dimension depends on various factors as discussed below. As the studied central point is by definition located in the middle of the testing cube, an isolated point, surrounded by values lower than the base of the cube s x s x s, generates a vertical line and a line generates a vertical plane. The slopes are respectively expressed by the following relations: y = -x, y = -2x with R² = 1. In the same way, a horizontal surface fills the half of the testing cube and the corresponding FD = 3 with R² = 1 (fig. 2). When the volume presents irregularities its fractal dimension decreases.

15The volume inside the cube depends on the coefficient h. The transformation into voxels can smooth the surface or, on the contrary, can enhance the surface roughness. When h = 1, the volume depends directly on the resolution of the DEM. Varying the coefficient h in Eq. (4) implies a change in the number of voxels vs, and thus the volume occupied inside the cube (fig. 1). High h values accentuate the surface roughness because the weight of the first term in this formula is more important than the second one (s/2). On the contrary, low values produce smooth surfaces because the second term of the formula is predominant. It is then possible to adapt the h value according to the nature of the DEM surface and the level of information expected.

16Determining the most appropriate range of grid sizes is a common problem for fractal dimension estimation (Foroutan-pour et al., 1999). The size of the frame used for the computation can be an even or an odd number; the range of grid size can be a power of two (Biswas et al., 1998) or a succession of integers. In our treatments, based on various tests applied to training zones, the size q of the sliding cube is chosen to be 12 or 24. Using these sizes, the slope P=ln(q)/ln(Ns) is more regular when exact dividers are employed. For a flat surface and a volume, the slopes P obtained are -2 and -3 respectively. Even if these proportions are respected and do not produce a different result after normalisation, odd sliding windows provide a slight deviation of the slope P for theoretical forms. The exact dividers, which define the range of grid sizes q, never correspond to half of the window side and then only a portion of the slope is taken into account in calculating the fractal dimension. When applying a 12 x 12 x 12 cubic pattern, five exact dividers (1, 2, 3, 4, and 6) can be found, excluding any border effect. In the second case (24 x 24 x 24), one can find seven exact dividers (1, 2, 3, 4, 6, 8, and 12). Treatments that employ these cubic sizes show that using the dividers 1, 2, 3, and 6 in the first case and the dividers 1, 2, 3, 6, and 12 in the second produces no bias. As the second term in Eq. 4 depends on the window size s, the use of the smaller window displays small irregularities; the larger window exhibits the general features of the DEM studied.

17The fractal dimension is calculated locally for each surface pixel in a sliding window centred on this pixel. As discussed above, the quality of results is greater with an even size than with an odd one. However, in the first case, the concept of ‘centre’ must be defined because the studied point cannot be located exactly in the centre but in one of the four points surrounding this centre. Using one of these points leads to the deviation of the result that favours one orientation.

18These remarks imply that the procedure must either calculate the local fractal dimension on each of these four pixels with the risk of blurring the response when calculating the mean, or define the value of the centre by taking the average value of the four pixels located at the centre of the window. Ic in Eq. (4) is calculated according to the second alternative in order to avoid blurring and to minimize the computation cost.

Test areas and results

19In order to test its ability to extract morphologic features according to different lithological and tectonic conditions, the local fractal dimension was applied in two regions. The first example, located in the region of Vittel (NE France), characterizes a sedimentary area affected by faults. The second example, located in eastern Anatolia (Turkey), is a volcanic region. In order to illustrate the performance of the procedure, the treatments concerning the Vittel region show mainly the variation of coefficient h and range of grid sizes effect, whereas the variation of the window size has been employed in the second example.

The Vittel region

20The Vittel area is located in northeastern France, 105 km west of the Rhine graben, between 48° and 49° N and 5° and 6° E. In the Vosges massif, outcrops of homoclinal Triassic and Liassic units overlie Variscan metamorphic and intrusive rocks. The east-west-trending Vittel fault is a major composite fault zone of complex geometry, resulting from the reactivation of a Palaeozoic discontinuity. In the region located between this fault and the N45-trending Esley fault, different groups of faults favour the water circulation that supplies the commercial Vittel spring water (Sykioti, 1994).

21The DEM of this area is an IGN (Institut Géographique National) product with a horizontal resolution of 50 m and a conical projection (fig. 3A). The corresponding geological map (Fig. 3B) has been geometrically corrected according to this projection in order to be overlaid on the DEM and compared with the local fractal results. The zone size is about 13.75 x 23.55 km (275 lines x 471 columns). Several treatments were applied to the DEM using a cubic size of 12 pixels and modifying only the value of the coefficient h. According to Eq. (4), decreasing coefficients contribute to smooth the roughness because a voxel inside the cube corresponds to a bigger hypsometric interval (ps/h). Then, the procedure detects only high altitude variations. On the contrary, when coefficient h increases, the surface roughness is accentuated and therefore small altitude variations are detected as well as large altitude ones.

22Thus, the position and extent of the Vittel fault (fig. 3C) are precisely highlighted and its prolongation westwards remains clear. When using high h coefficients (i.e., 100), every elevation difference is taken into account. Then, a thalweg feature in the testing cube corresponds to a filled cube with a gully whatever the altitude variation. The fractal dimension obtained is close to 3. On the other hand, a crest line is represented by a “wall” as formerly described (see Fig. 2) and its fractal dimension is close to 2 (fig. 3C). These results are obtained using a window size s = 12 and a range of grid sizes varying between 1 to s. These ensure the continuity of the features by producing a standardization of the estimated fractal dimension and increase the detection of crest and thalweg features.

23In contrast, the major units encountered in the studied zone were underlined using a low h coefficient but greater than one, because the region is relatively flat. With a coefficient equal to 2 (fig. 3D), the main structural and geological units are detected. A high value is obtained when, at the observation scale, the studied zone inside the moving window is relatively flat. The escarpment that corresponds to the Rhetian border is strong enough to induce in the moving window high roughness values codified with low grey tone values. Thus, the limit of the Rhetian appears clearly and corresponds exactly to the limits reported in the geological map. Concerning the Muschelkalk (limestones about 50 m in thickness) and the Bundsandstein (Vosgian sandstones), the incision depth of the drainage network is closely related to the respective nature of these two formations. For this reason, these units are clearly individualized and the contact corresponding to the N45-trending Esley fault is obvious, as well as the stratigraphical contact located at the base of Muschelkalk escarpment. The results are improved by superposing the geological map and the image resulting from the local fractal treatments. The frequency of the local fractal dimension depends on the sliding window size (fig. 4). Using a great size, the histogram is unimodal and becomes multi-modal with decreasing sizes, allowing us to obtain various classes corresponding to each mode.

The Ararat volcano

24With an elevation of 5123 m a.s.l., Mount Ararat is the largest volcanic centre and the highest point of eastern Anatolia. It corresponds to a compound stratovolcano formed by the ‘Greater’ and the ‘Lesser’ Ararat. This N150-trending, elongated massif lies inside a pull-apart basin with a similar trend (Pearce et al., 1990). This fault system is a horsetail splay structure and in addition to the right-lateral slip component, these faults have a normal component. The horsetail splay fault system cuts through the Greater and Lesser Ararat volcanoes and controls the position of the main eruptive centres of Ararat, as well as the alignment of parasitic cones (Karakharian et al., 2002). The volcanic activity of this region continued without interruption until historical times, possibly reaching its climax during the late Miocene and Pliocene (6 to 3 Ma). During Quaternary times, the volcanism appears to have been restricted to a few localised centres (Yilmaz et al., 1998). The youngest eruptions from the Ararat volcano are older than 10 000 years, and isotopic dating of the youngest lavas yields an age of 20 000 years (Yilmaz et al., 1998). Recently, geomorphic criteria applied to Mount Ararat have been used in order to study the tectonics of eastern Anatolia (Adiyaman et al., 2003). These criteria refer to the morphometric method developed by F. Garcia-Zuniga and J.-F. Parrot (1998) in order to define the base line and measure its elongation, and to the “Voxel wall” procedure (Baudemont and Parrot, 2000). The latter consists in parameterizing the DEM surface in a real 3D space using a vertical wall of voxels.

25The DEM of the Ararat zone (fig. 5A) is generated using the method developed by H. Taud et al. (1999). The pixel size equals 30 m. The size of this area is about 9.57 x 9.42 km (319 lines x 314 columns). The surface roughness of the different lava flows has been studied by means of different window sizes. A small window size (s = 12 with h = 4) reveals that the roughest surfaces mainly correspond to the more recent features, i.e., to the extent of the lava flows erupted from the recent parasitic cones located on the western flank. The extent of the eastern gully coming from the Greater Ararat is clearly observed, as well as the total base line of the whole edifice recently described by O. Adiyaman et al. (2003). A greater window size (s = 24 with h = 4) takes into account the regional structural features. As the local fractal values decrease in relation to brutal hypsometric variation as well as the character of the studied shape, a great window size reveals for instance the presence of isolated volcanic cones (low values), the presence of domes (high values surrounded by a crown of low values). The local fractal results depend on the size of the moving window as the frequency of the 2D signal is related to the surface roughness characterizing the type of studied lava flow. On the other hand, the value of the general roughness wavelength is an indicator of the regional tectonics. As illustrated by figures 5C and 5D, the SW quarter of the stratovolcano presents a globally lower roughness signature, the eastern limit of which emphasizes the presence of a large fault zone passing through the Greater Ararat crater. This fault corresponds to a branch of the horsetail splay structure described by A. Karakharian et al. (2002). In addition to the right-lateral slip, these authors assume that these faults have a normal component. The SW zone that has a weaker roughness response could indicate a huge collapse structure. In fact, the volcanic material (hypersthene andesite) observed in this zone (fig. 5B) is globally different from the material outcropping in the eastern part of the stratovolcano (mainly basalts and hypersthene basalts).

26Moreover, in order to study the effect induced by different window sizes and different values of h, the units of the geological sketch map drawn by Y. Yilmaz et al. (1998) have been reordered according to their petrographic nature. Six volcanic items (types of lava) are present: andesite, hyaloandesite, hypersthene andesite, basalt, hyalobasalts and hypersthene basalt. On the other hand, alluvium, fluvial and glacial fans, moraines and basement have been also taken into account in the statistical study. This geological map has been geometrically corrected in order to be overlain on the DEM and compared with the local fractal results. The mean of the local fractal dimension measured on the surface of each item is reported in the figure 6. The value of the fractal dimension depends on these two coefficients. Some items present specific values that allow us to characterize or to group them in different classes. When the window size is equal to 6 (fig. 6A), all the items except the alluvial formations are comprised in an unique class while the h value is equal to 1. For h values greater than 1, three classes can be observed: low fractal values are related to alluvium, basement and fans; median values to basalt, hyalobasalt and hypersthene andesite (the two formers are strongly connected); and high values to the remaining items. The same distribution is observed using window sizes of 24 where different values of h do not generate any variation (fig. 6C). On the contrary, with a 12 window size (fig. 6B) the items are classified differently. For instance, andesites are characterized by low values when s = 12 and high values when s equals 6 or 24. Therefore, using the mean value, statistical classifications can be obtained by taking into account the results produced by different values of coefficients s and h.

27The items corresponding to large zones seem to be formed by a heterogeneous ensemble corresponding to different erupted materials and the presence of different volcanic flows. High values are registered at the bulged top of the lava flows because this top appears as a flat surface at the observed scale, while the lateral flanks present lower values. This is the case for the isolated flows such as the hypersthene basalt from the Lower Ararat (fig. 5C). Hence, the high roughness alignments detected by the treatment, which are not present in the different geological units mapped by Y. Yilmaz et al. (1998), could be related to different volcanic flows running from the eruptive centres.

Conclusion

28This paper aims at examining the relation between the surface roughness and the geological and geomorphologic features. It is well known that the fractal dimension can be used as a parameter characterizing the surface roughness and the landscape shape. The treatment proposed here measures locally this dimension by an adaptive box counting 3-dimensional method. The algorithm is based on the use of two parameters: the size s of the moving window and a scaling factor h defining the hypsometric interval taken into account. This treatment has been applied to a DEM surface of two regions characterized by a difference concerning their lithological and structural conditions. The first example is related to the study of a sedimentary and faulted region in the Vittel area. In this case, two main results have been obtained according to the different configurations of the parameters involved in the calculation. As a first result, the limits of the geological units detected by the local fractal dimension correspond to these limits as they have been mapped on the field. The position of these limits match either to the border of a homogeneous morphologic feature or to the transition zone between two different incision depths of the drainage network. On the other hand, by using low hypsometric intervals by means of the coefficient h, the procedure can detect fine features. We can therefore underline the Vittel fault as well as its western prolongation. The second example focuses on the volcanic region of Mount Ararat. In this case, the relation between the local fractal measurements and the different volcanic formations which have built up this stratovolcano has been studied by using a statistical approach. Different volcanic classes can be distinguished by using their mean fractal dimension value. On the other hand, the results provided directly by the local fractal dimension show the extent of the volcanic flows according to their surface roughness. Furthermore, a collapse structure related to the regional strike slip faulting has been detected. The results obtained with the two training DEMs corresponding to lithological conditions that are completely different, corroborate the efficiency of the procedure. The latter allows the retrieval of information about small structural features in flat zones as well as general morphologic characteristics of rugged areas.

29Based on the first results, one can conclude that the roughness of a surface is strongly correlated with the nature of the material that forms the geology of a studied region. In the case of the Vittel region mainly formed by a sedimentary sequence and studied as a training set, the application of the treatment by means of a descriptive analysis shows that this technique is efficient to extract and precise the meaning of the different morphometric features. In the case of more complex region such as the Mount Ararat that presents unknown structures, the extraction and the meaning of the morphometric features resulting from the treatment, needs to understand the nature of the filtering result; with such a goal, it is necessary to analyse precisely the results provided by using different values of the coefficients s and h, in order to realize an objective synthesis.

30Actually, the local fractal dimension provides useful information about geological and geomorphologic features. The proposed method detects different types of structures according to the observation scale and provides useful information concerning the geological and geomorphologic mapping.