Spatial Variation of Fractal Parameters and Its Geological Implications

The spatial variation of the fractal dimension and its geologic implica­ tions are studied using an empirical approach. The fractal parameters, including the surface's fractal dimension, D .. the angular fractal dimensur[� sion, D , the gamma ( y) value and the break-distance (R), are derived by ang the variogram method. Synthetic surfaces are generated by the successive random addition (SRA) method by assigning different values to the Hurst exponent. Two test sites were selected to study the link between the fractal parameters and the geologic features. The following conclusions can be drawn from this research: (1) The fractal parameters of a landscape sur­ face are region-dependent and scale-dependent. The spatial distribution of fractal parameters should be studied in a suitable size of spatial unit. (2) The angled variogram method can readily disclose the anisotropic nature of a landscape surface. (3) The surface's fractal dimension reflects the litho­ logic variations underlying a landscape surface, and the gamma value re­ flects the topographic relief of the surface. (4) The mean direction and the vector resultant of angular fractal dimensions have a close relation with the major geological structures of the landscape. (


INTRODUCTION
The use of fractal geometry and, in particular, a fractal dimensions (D) to describe the nature of landscapes is now a growing field of research. The theory is based primarily on the work of Mandelbrot (Mandelbrot, 1982;Mandelbrot, 1975). Many parameters have been developed for characterizing landscapes and related properties. The fractal dimension of a surface appears to contain useful information about the surface which is not given by other morphometric measures (Feder, 1988;Klinkenberg, 1992). The fractal model treats the land scape region being investigated as a distinct entity and the fractal nature of the region is ini tially assumed to be homogeneous (Xu et al., 1993). Therefore, one can use mathematical functions to describe certain characteristics of the landscape region and should be able to 1Department of Earth Sciences, National Cheng Kung University, Tainan, Taiwan ROG obtain a uniform fractal dimension. However, landscapes are a composite of many competing and complicated geological processes, such as faulting, folding, flexure, erosion and sedimen tation. Milne (1991) indicated that exact fractal patterns are unlikely to occur in landscapes because contemporary patterns are the result of several processes that dominated in the past. Therefore, a landscape should be called statistically self-similar or statistically self-affine if it possesses a fractal nature (Xu, et al., 1993). Research to date indicates that self-similarity is exhibited only in limited regions and over limited ranges of scale in real landscapes. Most studies more or less confirm the spatial variation of D over landscape surfaces, namely the scale-dependent variations and the region-dependent variations (Burrough, 1981;Mandelbrot, 1977Mandelbrot, , 1982. But, as yet there is no explicit criterion for defining how large the region with homogeneous fractal properties for which the landscape is statistically self-similar needs to be, or for how to delimit the region (Xu, et al., 1993). Indeed, some doubts remain about the estimation of D with small samples other than the currently used 256 X 256 data matrix.
There are many methods for determining the fractal dimensions of landscapes. The fractal dimensions computed using different methods for the same region of a landscape are often different (Xu et al., 1993;Gallent et al., 1994;Klinkenberg and Goodchild, 1992). Having so many methods of computing the fractal dimension is inadequate because it could conceal problems, such as differences in their mathematical basis and mechanisms as well as under standing of the physical significance of fractals (Xu, et al., 1993). The links between the fractal dimensions and the physical processes which produce that characteristic of form cap tured by D have not been identified as of yet (Klinkenberg, 1992). This is an important issue which remains to be tackled before the concept of fractal dimension can be widely accepted among experts in the geologic society. Nevertheless, Klinkenberg and Goodchild (1992) pointed out that the variogram methods produced parameters which appear to discriminate between physiographic provinces, and that even within physiographic provinces some differentiation is possible. This implies that the methods produce robust and consistent results. Furthermore, the oriented (angled) variogram methods may extract the anisotropy of the terrain feature that reflects the regional structures of valleys and ridges. Although the accuracy of the digital elevation model could be an influential factor in scale-dependent variation and directional bias of fractal dimension, the angled variograms are picking up detail that is not visible in a contour plot (Klinlenberg and Goodchild, 1992).
Taiwan is situated in the orogenic belt of the Cenozoic arc-continent collision between the Eurasian continental margin and the Luzon volcanic arc of the Philippine Sea plate (Ho, 1982;Teng, 1990). The major stratigraphic units and the major geologic structures are distrib uted along the elongated arc-shaped framework trending north-southerly. Topographic fea tures are highly controlled by this geologic framework and exhibit an apparently anisotropic nature in the regional scale. It is relevant to study the spatial variation of fractal parameters using the digital elevation models (the DEMs) of landscapes on Taiwan. The objectives of this research are as follows: (1) We try to define the smallest spatial unit within which stable fractal parameters can be computed by the variogram method. (2) Synthetic surfaces are generated by the successive random addition method (Voss, 1985) in order to testify if the spatial variation and the directional bias of the fractal dimension are caused by the sample size used in the estimates of fractal dimension. (3) After the proper size of the DEM data matrix has been empirically determined, a couple of test sites, which have been geologically mapped, are used to explore the links between the fractal parameters and the geological features. The geological implications of the fractal parameters then are fully discussed.

FRACTIONAL BROWNIAN SURFACE
Until now, Mandelbrot's fractional Brownian motion (fBm) or fractional Brownian func tions probably provide the most useful mathematical model for the random fractals found in nature, particularly for application to landscape analysis (Xu, et al., 1993). In practice, frac tional Brownian motion is a random process F(t) with Gaussian increments and where H, the Hurst exponent, varies from 0 to 1. When properly rescaled, the two random functions F(t) and F (rt) are statistically similar. For a surface, the single variable ti s re-rH placed by point coordinates x and yo n a plane to give F (x, y) as the surface altitude z at position x, y. The surface that consists of these F (x, y) points is usually called a fractional Brownian surface. On a fractional Brownian surface, the variogram is described by The variogram takes on the form of a power function in which H should vary between 0 and 1.
In the case of profiles, the fractal dimension, D, has the following relation with H, D=2-H while for fractional Brownian surface,

D=3-H
As H increases toward its upper limit (i.e., small D), the variability of the surface is small locally, but rises rapidly with distance, whereas, when His small (i.e., large D) the surface shows high local variability but a slow increase at large distances.

SIMULATION OF FRACTIONAL BROWNIAN SURFACE
Several algorithms have been proposed to simulate a fractional Brownian motion and a fractional Brownian surface. A summary of the algorithms is given in Voss (1985). We have used a successive random addition (SRA) method (Voss, 1985). The SRA method generates fractional Brownian surface having a power-law variogram model. In this experiment we modified the SRA method by introducing a different value of parameter H while adding the midpoints of the square lattices. In the mesh pze, 8 denotes the resolution of such a grid. We obtain another square grid of resolution o 10 2 by adding the midpoints of all squares. In this case, the orientation of the new square lattice is rotated by 45 degrees. Again adding the midpoints of all squares gives us the next lattice with the same orientation as the first one, :µid the resolution is now 8 12. In each stage, we thus scale the resolution with a factor of r =11 0 2 , and we add random displacements using a variance which is r 2 H times the variance of the previous stage. Different values of the parameter, H, can be introduced into the 90 degrees' lattice and the 45 degrees' lattice respectively. A fractional Brownian surface with anisotropic nature (i.e., with directional fractal dimension) can be generated in this way.

ESTIMATION OF FRACTAL DIMENSIONS FROM VARIOGRAMS
The variogram method is a well-known technique for determining the spatial characteris tics of data (Kridge, 1966;Agterberg, 1982). The essence of this technique is that the statisti cal variation of the elevations between Salnples is some function of the distance between them.
The independent variable is the distance between pairs of points, the dependent variable is the (semi-)variance of the differences in the data values for all samples the given distance apart.
By making the assumption that land surfaces have statistical properties similar to those of fractional Brownian surface, it is possible to obtain the fractal dimension of the landscape surface directly from the variogram (Mandelbrot, 1975;Mandelbrot and Van Ness, 1968).
The fractal dimension is directly related to the slope of the best-fitting line produced when the log of the distance between samples is regressed against the log of the mean-squared differ ence in the elevation for that distance. The 95%H confidence interval of the best-fitting line is calculated to estimate the accuracy of the fractal dimension. The intercept of the best-fitting line with the ordinate is define<l as the gamma ( y ). The break-distance (R) denotes the dis tance within which the statistical �elf-similarity can be applied. Fractal dimension, D, y, and Rare the fractal parameters estimated from the variogram method (Figure l).
In this research, the angled variogram method was applied to the DEM subsamples with spatial resolution of 40m. The angular fractal dimension (D ) was calculated along seven ang azimuths using the following sampling schemes: (1) Elevation pairs were constrained to fall in the same row or column of the DEM subsamples; the associated dimensions will be referred to as D0 and D90 to represent the north-southerly and the east-westerly directional fractal dimension.
(2) Elevation pairs were constrained to fall in the diagonal of the DEM subsamples; the associ ated dimension will be referred to as D45 and D135 to represent the north-easterly and the south-easterly directional fractal dimension.
(3) Elevation pairs were constrained to fall in the direction which is aligned by Z.. and Z. 2 , . 1, IJ 1+ J· or by Zi i and Z;+ 2 , j+t; the associated dimensions will be referred to as D26 and D154 to represent the north-northeasterly and the south-southeasterly directional fractal dimension.
(4)Elevation pairs were constrained to fall in the direction which is aligned by Z. and Z. 1,. 2, or IJ I+ J· by zij and zi+ I ' j+2; the associated dimensions will be referred to as D64 and Dll6 to represent the east-northeasterly and the east-southeasterly directional fractal dimension.
(5)The entire elevation pairs were sampled to calculate the entire surface dimension; the asso-0.01 dated fractal dimension will be referred to as D surf.

659
The fractal rose diagram can be constructed using the D as the radius and the angle as ang the azimuth. The mean direction and the vector resultant are also computed to demonstrate the anisotropic nature of landscape surface. The directional fractal dimension is deliberately ex pressed by the vector format with azimuth ranging from 0° to 180°.

RESULTS AND DISCUSSIONS
A series of land surfaces, each having 1024 by 1024 arrays, were synthesized by the SRA method using the Hurst exponent, H=0.2, 0.5 and 0.8 respectively. Since the synthesized land surfaces should exhibit behavior consistent with fractional Brownian models, all variograms should produce similar dimensions. The synthesized land surface should be isotropic in na ture, so D rr should be uniform even when the land surface is dissected into smaller units. In SU addition, a series of land surfaces, with the same size, were also synthesized by the SRA method using different values for H" . in this way, should be anisotropic in nature and a directional bias of the fractal dimension should be expected. However, note that the SRA method may not be an appropriate method to generate a fractional Brownian surface, arid the variogram method may not be an accurate estimate to the surfaces' fractal dimension (Wen, et al., 1997). Nevertheless, we decided to use these methods because they are relatively easy to implement in our computer system, and because we are mainly concerned with the geological implications contained in the fractal parameters, not in comparing different estimating methods. The following issues are the ma jor concerns in this research and will be discussed in detail: (1) From how small a data matrix can the fractal dimension be reasonably estimated?
As noted by Wen and Sinding-Larsen (1997), the variogram method results in overesti mates .for small values of the underlying fractal dimension and underestimates for large values of the fractal dimension. Furthermore, the precision of the estimates is dependent both on the number of samples used in the fitting of the variogram and the underlying "true" value of the fractal dimension to be estimated. Using variogram values calculated from large lags (i.e., larger data matrices) leads to a reduction in both the accuracy and the precision of fractal dimension. But, how small a data matrix is appropriate for a reasonable estimate of the underlying fractal dimension? We did the following experiments to determine the proper size of the data matrix. We first generated three sets of fractional Brownian surfaces with 1024 by 1024 data matrix, each of which had the theoretical fractal dimension of 2.2, 2.5 and 2.8 respectively. The fractal dimension of the synthesized surface:-.. was calculated by moving a window over the entire surface using the variogram method. The window size and the offset for each move could be selected arbitrarily.
We defined that the calculated fractal dimension which fell within the range of ±0.l of the theoretical fractal dimension was an acceptable estimate, because a variation of ±0.1 was good enough to tell the difference of the fractal dimension among 2.2, 2.) and 2.8. The plot of window size versus percentage of acceptance showed that the percentage of acceptable esti mates decreases gradually as the wind ow size becomes smaller and it drops drastically as the window size is smaller than 30 on each side (Figure 2). We suggest that 30 by 30 be the smallest data matrix from which one can obtain about 80% of accuracy in the estimate of the surfaces' fractal dimension. Referring to the conventional DEM available in Taiwan, the smallest spatial unit is l .2km x 1.2km for estimating the fractal parameters.
(2) Does the angled variogram method produce directional bias in estimation of fractal dimen sion?
Few studies have considered the anisotropic nature of the fractal dimension of landscape surfaces directly. Klinkenberg and Goodchild (1992) found that certain physiographic prov inces, the Basin and Range, for example, are much more heterogeneous in their angular fractal characteristics than others. However, Burrough (1981) had pointed out that the form of the variogram is often highly dependent on sampling direction and sampling interval. When a sampling interval matches a particular scale of a phenomenon in the landscape, it perceives an apparently lower fractal dimension.
It is, therefore, important to clarify if the angled variogram method results in the direc- ( c) The SRA method may not have generated a random fractional Brownian surface; neverthe less, the angled variogram method still can disclose the anisotropic nature of the synthetic surfaces in a relative sense ( Figure 6).
Therefore, it may be concluded that the angled variogram method applied in a moving window fashion can be used in the study of the spatial variation of the fractal dimension over a landscape surface.  (3) How to interpret the fractal parameters in the geological sense?

663
The fractal parameters, which can be derived by the variogram methods, include the fractal dimensions (D,u,r and Dang), the ordinate -intercept ( y ), and the break -distance (R). The vector resultant and its mean direction can be calculated from each angle fractal dimension.
Various authors have attempted to relate differences observed in the fractal dimension of vari ous surfaces, and within various surfaces, to a number of geomorphic factors. (Ahnert, 1984;Culling, 1986;Culling and Datko, 1987;Fox and Hayes, 1985;Chase, 1992;Klinkenberg and Goodchild, 1992). It is too sweeping a statement to say that the fractal dimension character izes only "irregularity" or "roughness" of the terrain (Xu et al., 1993). A landscape may be the result of a number of opposing processes, with each process resulting in a specific fractal domain (Klinkenberg and Goodchild, 1992). On the other hand, some experimental studies have shown that the gamma value has a distinctively positive relation with the relief index of a topographic surface (Chao, 1995).
In this experiment, two test sites were selected to study the links between the spatial variations of fractal dimension over landscape surfaces and the geological factors. They are located at the southwestern foothills of Taiwan. Imbricate folding and thrusting of Quaternary The angled variogram method was applied to the DEMs by the moving window operation.
The window size of 33 rows by 33 columns with an offset of 16 rows by 16 columns was used to acquire a higher resolution in the fractal dimension to discuss the spatial variation of fractal dimension over the test sites. The following are some observations made in the experiment: (a) Test site A (Figure 7), sedimentary rocks in the west. In addition, several west-vergent major folds, whose axial planes also trend northeasterly, are bounded by the major thrust faults in the imbricate system.
The contour map of the surface fractal dimension (Dsurr) exhibits a pattern controlled by the lithology of the area (Figure 7). The Dsurr of the sandy sequences has a lower value of 2.26±0.08, while the D r of the shaly sequences has a higher value of 2.34±0.1. The D rf of w the alternated sandstone and shale has the value between the previous two. The Dsurr of the submetamorphic sequences has the lowest value of all (Table 1). The pattern of Dsurf contours strikes northeasterly and coincides with the ridge and valley structure. The contour map of the gamma value also exhibits a striking coincidence with the ridge and valley structure. The gamma value of the ridge system is relatively higher than that of the valley system. In short, the ridges, where the thickly bedded sandstone crops out, are characterized by a low D su rf and a high r , while the valleys are vice versa.
The spatial distribution of the vector resultant and the mean direction of the angular fractal dimension exhibit a distinct pattern controlled by the local structures of the area (Figure 8). Generally speaking, the mean direction of the angular fractal dimension follows closely the major structural trend in this area. The length of the vector resultant indicates the degree of anisotropy of landscapes. We can observe the landscapes, where constructed by the massive sandstone, exhibit more anisotropic characteristics than those constructed by the shale. More over, the vector resultant and the mean direction of angular fractal dimension change dis tinctly across the major faults and folds. This implies that topographic features in the Ghia -Hsien area are highly controlled by the geological structures and can be readily depicted by the fractal parameters.
(b) Test site B (Figure 9), the Chia-Yi area, is an area which shows a drastic change in lithology as well as in topography, when contrasted with test site A. A major thrust, the Chu Kou fault, separates the Neogene elastic fades of the imbricate fold and thrust belt in the east from the Quaternary mollasse fades in the west. It also marks a major physiographic bound ary between the rugged mountainous terrain, and the rolling hills and the coastal plains. The major thrusts trend northeasterly and are cut by a few northwest and eastwest trending transcurrent faults. Some of the faults are probably still active. The Neogene elastic fades  14. The D,0rr of the coastal plain has a surprisingly high value of 2.67±0.14 ( Table 2 ). A contour map of the gamma value of the test site B also matches the topographic features ( Figure 10). The gamma value reflects not only topographic elevations, but also topo graphic relief. Spatial distribution of the mean direction and the vector resultant of the angular fractal dimension generally follows the major structural trend of the area except for some local variations( Figure 8). The local pattern, strongly influenced by the geological structures, im plies that these structures may occur relatively young in age and can be readily disclosed by the angled variogram method.  The following conclusions can be drawn from the above experiments: (1) The surface's fractal dimension, Dsurr' and the gamma value both are useful parameters for differentiating rock units in different geological provinces. For example, the coastal plains are characterized by a high D,urrand a low y, while the sedimentary rocks in the imbricate fold and thrust belt are characterized with a low D,u,r and a high y. The realization of the contrastive features simulated by the SRA method is shown in Figure 11. This seems to agree with the observations made by Pentland(1984), and by Klinkenberg and Goodchild(l 992) that the surface's fractal dimension reflects the roughness of landscape surface and the gamm a reflects the topographic relief However, we propose that the smface' s fractal dimension be a measure of texture of a landscape surface. "Texture" is usually defined as frequency of change and arrangement of tones in the field of image processing and photo interpretation. In the geomorphologic sense, texture can be defined as the fre quency of change and arrangement of topographic heights. It means something more than what roughness can tell about a surface. High fractal dimension implies a quick change of elevations in a local area but a slow change of elevations at a large distance, while low fractal dimension implies a small local variation in elevation but a large variation at a long distance. This can well explain the fractal characteristics of the coastal plain and the moun-  Fig.11. The synthetic surfaces generated by the SRA method. (a) The realiza tion of the Coastal plain, generated using H=0.2 (D=2.8) and relief=50m, has a low, flat but a rough surface. (b) The realization of the mountain range, generated using H=0.8 (D=2.2) and relief=800m, has a higher relief but a smoother surface.
tain ranges in the test sites.
(2) The angular fractal dimension can readily depict the anisotropic nature of landscape sur face, which is controlled by the geological structures. Therefore, the geological structures can be disclosed by the angled variogram method. Which geological processes cause the directional bias of fractal dimension is to be explored. For the case of the Ghia-Hsien area, the mean direction of the angular fractal dimensions is apparently parallel to the major thrusts and folds. It may suggest that the extensional jointing developed perpendicularly to these thrusts and folds cause a high frequency of change in elevations (i.e., high fractal dimension).
(3) The break-distance seems to have little relation with topographic features and geological structures. It is constrained by the size of a moving window. Its geological implication has not been determined by this research.

CONCLUSIONS
The spatial distribution of fractal parameters over a landscape surface and its geological implications have been studied by an empirical approach. Several experiments were done on the synthetic surfaces generated by the SRA method, and the angled variogram method was used to calculate the fractal parameters of a landscape surface. The following significant re sults were obtained by this study: (1) A landscape surface is anisotropic under geological control. The fractal parameters of a landscape surface are region-dependent and scale-dependent. The spatial distribution of fractal parameters can be studied in a moving window fashion. Accuracy and precision in the estimation of fractal parameters are influenced by the size of spatial unit.
(2) The angled variogram method can readily disclose the anisotropic nature of a landscape surface.
(3) The surface's fractal dimension reflects the lithologic variations underlying a landscape surface, and the gamma value reflects topographic relief of the surface.
( 4) The mean direction and the vector resultant of angular fractal dimension have a close relation with the major geological structures.