## Abstract

Surface roughness is an important rock property that is measured for structural geology and engineering purposes. We have developed an automatic technique to detect anisotropic features on rock faces based on fractal analysis. The analysis method has been applied to synthetic surfaces, and to digitally mapped point clouds of natural rock surfaces shaped by weathering, fault wear, and mining. We illustrate the technique using field examples from Permian sandstones containing brittle shear zones in northeast England, the surface of a neotectonic fault in Turkey, Proterozoic quartzite from central Australia, and Devonian Quartzite in an aggregate quarry in Germany. Roughness analysis of these natural examples suggests that a significant change of roughness value, anisotropy, and anisotropy direction can exist across scale. Our analysis method represents a step toward developing a toolkit to automatically detect and interpret surface characteristics from digitally acquired data sets. It has widespread potential for applications in rock engineering and the geosciences.

## INTRODUCTION

Geometrical characteristics of rock surfaces are an important source of information in structural geology and engineering. One of these surface characteristics is roughness, which is a surface's deviation from planarity at a defined scale of observation. Rock surface roughness has been described as a self-affine fractal feature because its statistical properties are retained on different scales in different directions (e.g., Brown and Scholz, 1985; Kulatilake and Um, 1999). Because the profiles of many natural rock surfaces are approximately self-similar within small wavelength bands, fractal methods are often used to measure a roughness along profiles (e.g., Power and Tullis, 1991).

In geotechnical engineering, surface roughness of joints is an important material property because it influences the failure behavior of a fractured rock mass. For geotechnical purposes, the scale of roughness is of importance. On a small scale (millimeters to centimeters), roughness on fractures and joints will determine shear strength, whereas on a larger scale, roughness, in this case sometimes referred to as waviness, will influence the direction of shear and the dilation of the rock mass during failure [e.g., International Society for Rock Mechanics Commission on Standardization of Laboratory and Field Tests (ISRM), 1978; Hobbs, 1993; Grasselli et al., 2002].

Deformation structures such as faults, shear zones, fractures (including veins and joints), foliations, and lineations are often expressed as morphology on rock surfaces, and their shape and orientation are an important source of information in structural geology (e.g., Passchier and Trouw, 2005). Rock surface roughness is often anisotropic, because planar discontinuity sets intersect the surface, or, in the case of exposed fault surfaces, structural anisotropies have developed due to shearing on the surface (Anderson, 1951). Geometrically irregular fault surfaces are of great significance to understanding locally inhomogeneous stress states along active faults (Chester and Chester, 2000). A number of studies have used manual measurements to show that the roughness of fault surfaces is greatest in the direction perpendicular to slip, and smallest in the direction parallel to slip (Hancock and Barka, 1987; Power et al., 1987, 1988; Power and Tullis, 1991, 1992; Lee and Bruhn, 1996; Jackson and McKenzie, 1999). As McCaffrey et al. (2005) pointed out, digital mapping techniques can greatly facilitate data acquisition and processing for many geoscience purposes, because digitally acquired data can easily be stored, shared, imaged, and analyzed in a spatially referenced manner. Digital mapping of rock faces using light detection and ranging (Lidar, or laser scanning) and digital photogrammetry methods require only limited access to the rock face, and are much more time efficient than labor-intensive manual measurement procedures. Based on digitally acquired data sets, Sagy et al. (2007) showed that roughness parallel to the transport direction decreases as fault slip length increases, Jones et al. (2006) related heterogeneous slip patches of the Arkitsa fault in Greece to its surface geometry, and Renard et al. (2006) reported that striations on a fault surface in the French Alps displayed self-affine scale invariance.

In this paper we present a new workflow in using a fractal analysis method to measure and plot roughness along different orientations. We also show how this method can be applied to measure surface anisotropy on rock surfaces shaped by erosion, fault wear, and blasting. Our method allows the comparison of surface anisotropy across scales, and variations of roughness with direction. It also allows the isolation of rock surface features of a specific range of wavelengths with a careful selection of the scales analyzed. The proposed system is similar to those of Renard et al. (2006) and Thomas et al. (1999) in that it analyzes the anisotropy of a surface. However, Renard et al. (2006) used a statistical analysis method using a root mean square correlation measure for the one-dimensional (1D) and 2D analysis of the surface. Thomas et al. (1999) analyzed the surface anisotropy using a structure function based on autocorrelation between the surface points for their tribology research.

## ROUGHNESS MEASUREMENT WORKFLOW

Methods for measuring surface roughness differ depending on the purpose of measurements and the materials investigated. When studying friction, lubrication, and wear of engineered materials, engineering tribologists commonly use central line average roughness, root mean square roughness, mean value of the maximum peak-to-valley height, and ten-point average (Stachowiak and Batchelor, 1993). Similar methods are also used to investigate Earth surface processes and landform evolution (McCarrol and Nesje, 1995).

The roughness of a rock face can be measured by considering profiles in different orientations and positions (ISRM, 1978). These profiles are used to determine how much the roughness varies across the surface and how isotropic the surface is. Statistically, this means taking a sample population from a continuous population space, the rock surface. The roughness along any given orientation is represented by calculating the mean and standard deviation of roughness measurements for all profiles of the same orientation. Statistics from sample profiles taken at different orientations are used to measure the anisotropy of the surface along with calculating the directions of maximum and minimum roughness.

We apply our analysis to point cloud data generated from three-dimensional photogrammetric surface models. The software package used to generate these surface models is Sirovision, a commercial code that is mainly designed for geotechnical purposes (Poropat, 2001; www.sirovision.com). Sirovision generates three-dimensional surface models from a pair of digital images, which are either taken from surveyed camera positions, or contain surveyed points in the images. The 3D surface models consist of triangular surfaces derived from mul-tipixel facets. The point clouds analyzed in this study represent the vertices of these triangular surface elements.

The output of our analysis is a polar plot that represents surface roughness on varying orientations (Fig. 2). Given the point cloud data, a number (N) of profiles for each orientation is sampled from the data. Then a roughness measure for each of the N profiles is calculated, and the mean roughness value for the N profiles is computed and used to represent the roughness of the specified orientation. The specific steps in this analysis include:

Transform the points of the rock face so that it faces straight up and the largest width of the face is aligned with the

*x*axis using principal components analysis (PCA).For each number (R) of rotations, rotate the surface about the z axis to a specified orientation: sample N profiles that are parallel to each other but all perpendicular to the plane of best fit (Fig. 2B); calculate the roughness for each of the N profiles; and calculate the mean and variance of the N roughness measurements.

Calculate the anisotropy.

Identify the direction of maximum and minimum roughness.

Visualize the analysis result.

### Rotating the Surface

The rock surface is represented as a point cloud, on which we perform a PCA. The objective of the principle component analysis (PCA) is to line up the normal of the plane of best fit to the *z* axis, such that the direction of maximum variance of the point cloud representing the surface is aligned with the *x* axis, the direction of intermediate variance with the *y* axis, and the direction of minimum variance with the z axis. Figure 2A gives an example of how the new axes align with the original rock surface. In this diagram the longest cylinder will become the *x* axis. The steps followed in this analysis are:

Center the point cloud by subtracting the mean point from each point within the point cloud.

Calculate the point cloud's covariance matrix by multiplying the list of points by its transpose.

Calculate the ordered eigenvectors of the covariance matrix (the eigenvector associated with the largest eigenvalue is aligned with the direction of maximum variance).

Multiply the list of points by the eigenvalue matrix.

### Sampling of Profiles

The model surface is sampled along parallel profiles within the *y-z* plane, while the surface is rotated about the surface normal (the z axis) using a transformation matrix (Fig. 2B). The roughness of each profile is calculated and the measurements for all profiles are used to compute the mean and standard deviation of the roughness for the sample (orientation).

### Roughness Measurement

The fractal dimension of a profile is a known measure of roughness (ISRM, 1978). The fractal dimension of a curve is calculated by comparing the curve length when sampled at different scales (Barnsley, 1988). Theoretically the fractal dimension is a limit calculated across all scales, but realistically the scales are limited by the data. In our analysis we have used the compass walking technique (Bourke, 1998; Gemperline and Siller, 2002). In this method the length of the curve is measured by how many times the compass intercepts a curve across its length. Then the compass length is lowered by dividing its current length by 2. Reasonable values for the absolute minimum and maximum compass lengths depend on the size and the spatial resolution of the model: the minimum compass length is limited by the resolution, and the maximum compass length is limited by the size of the model. Realistically, the minimum compass length is further limited by the precision of the point cloud. The sensitivity of the calculation to the minimum and maximum compass lengths relative to the feature of interest can be seen in Figure 1. Here the fractal dimension of a sine wave has been computed with a compass length of 3.32–0.11 units as the wavelength and amplitude of sine wave vary.

### Visualizing the Result

The results are visualized in a polar plot. The theta value of each point in the plot corresponds to the angle the samples are taken, as measured in Figures 2B and 2C. The distance from the origin of each point corresponds to an adjusted roughness value for that orientation (Fig. 2C). The adjusted value is calculated by taking the mean roughness value of the orientation and subtracting 1. This can be done because a profile always has a fractal dimension between 1 and 2. The use of the adjusted roughness value helps to accentuate differences in roughness on the polar plot since values are always larger than one.

For the analysis purpose, we identify not only the minimum and maximum directions of roughness, but also the directions that show significantly large or small roughness measures. This is achieved by the following processes. The roughness values are plotted against the orientations, and then a central moving average is applied to smooth the data. The amount of smoothing is controlled by the number of points used for each average. We found that 12 points worked well for a sample of 120 orientations. This smoothed curve was then used to find maxima and minima extremes. The largest and smallest of these were used to specify the directions of maximum and minimum roughness.

The surface anisotropy is measured by calculating the roundness *a* of the rose plot, where *a* has range of 0 < *a* < 1, and 1 would be a circle representing perfectly isotropic surface roughness.

## EXPERIMENTAL RESULTS

### Roughness Analysis of Synthetic Surfaces

Three synthetic surfaces based on a corrugated iron geometry have been generated, each with different amplitudes and wavelengths. Note that these surfaces are not fractal surfaces, but are used to understand the effects of different compass lengths on our roughness analysis. The first surface has a wavelength of 1 unit, an amplitude of 0.15 units, and wave crests and valleys are oriented parallel to 0° (Fig. 3A). The second surface has a wavelength of 0.4 units and an amplitude of 0.05 units, but wave crests and valleys are oriented at 45° (Fig. 3D). The third surface is the interference pattern that results from combining the above surfaces.

These surfaces were analyzed to determine the direction of maximum roughness, *R*_{max}, and the roundness (i.e., the isotropy measure). When a long-wavelength–high-amplitude surface (Fig. 3A) is analyzed using a compass length range of 0.64–0.11, *R*_{max} is 0°, and the roundness is 0.26 (Fig. 3B). The direction of maximum roughness is correctly detected, and the roundness value demonstrates that the surface is very anisotropic. However, when the same surface is analyzed with a compass length range of 0.11–0.02 (only small compass lengths), the analysis detects *R*_{max} as 38.1° (Fig. 3C), which is incorrect. Measurement at this scale produces an unreliable detection of roughness orientation, because the compass range does not contain larger compass lengths that are necessary to measure the roughness of this surface. A short-wavelength–low-amplitude surface (Fig. 3D) is analyzed in a similar fashion. When analyzed with a compass length range of 0.64–0.11, *R*_{max} is 24°, and the roundness is 0.34 (Fig. 3E). For a compass length range of 0.11–0.02, *R*_{max} is 46°, and the roundness is 0.17 (Fig. 3F). In this case the direction of maximum roughness is approximately correct for both of the ranges as the smaller wavelength of the surface is detected by the small compass lengths. When the combined surface (Fig. 3G) is analyzed using small compass lengths, the statistics of the short-wavelength–low-amplitude surface are isolated (cf. Figs. 3F and 3I). Here *R*_{max} is 50° when using a compass length range of 0.11–0.02, which is within 5° of the correct direction. A compass length range of 0.64–0.11, however, results in the detection of *R*_{max} at 23° (Fig. 3H) with a broad spectrum of directions of high roughness values near 23°. This is due to the statistics of combining two surfaces, generating a wide range of directions where profiles are rough.

These results show that our methodology can accurately determine values of *R*_{max} and roughness isotropy, but that these values are strongly dependent on the choice of an appropriate compass length range in relation to the length scale of the surface feature.

### Roughness Analysis of Field Data

The rock surfaces we analyzed were derived from point clouds generated in Sirovision, (Poropat, 2001; www.sirovision.com). For our analysis, we have chosen four field examples with significantly different surface characteristics: (1) Permian sandstones with conjugate brittle shear zones in Cullercoats, northeast England; (2) the footwall surface of the Yavansu fault in Turkey; (3) Proterozoic quartzite with spaced fractures in central Australia; and (4) Devonian Quartzite in an aggregate quarry in Germany. The outcrops are described in more detail below, along with the results of the roughness analysis.

#### Permian Sandstone at Cullercoats

Permian sandstone outcrops at Cullercoats, county Tyne and Wear, UK, provide a good example for anisotropic surface roughness. At Cullercoats the 90-fathom fault displaces a hanging-wall sequence of mid-Permian dolo-stones and sandstones against a footwall of Carboniferous coal measures. In the sandstones, deformation related to the movement of the fault is expressed as a set of cataclastic brittle shear zones, also termed deformation bands. These deformation bands are well cemented, and therefore create a positive relief, which makes them a dominant surface feature in outcrops weathered by wave action. While the deformation bands appear to represent a bimodal distribution (conjugate set) in our model, the actual fracture geometry across the entire outcrop is interpreted to be of quadrimodal symmetry (De Paola et al., 2005; Healy et al., 2006).

Our surface model represents an ∼0.8 m × 0.6 m area of an approximately west-southwest–east-northeast–oriented vertical rock face at Cullercoats beach (55°2′1.43″N; 1°25′53.90″W). The large-scale analysis with a compass length range between 2 cm and 0.35 mm detects high roughness values at 30° and 156°, a high level of roughness around the image horizontal, and a relatively low anisotropy (Fig. 4C). This is in good accordance with the orientation of the brittle shear zones perpendicular to these orientations (Figs. 4A, 4B). At a compass length range between 3.5 mm and 0.6 mm (Fig. 4D), the anisotropy is dominated by smaller scale features. Likely candidates include sedimentary bedding, which is dipping steeply toward the left (Fig. 4A), and other short-wavelength linear features at the lower left and upper right of the image. Another difference to the larger compass length range measurement is that at the small scale there is much less roughness in the horizontal.

#### Yavansu Fault at Kuşadasi

The footwall surface of the neotectonic Yavansu fault displays surface features on different scales. The Yavansu fault is part of a post-Serravalian to recent set of east-west–trending normal faults that accommodates the north-south extension in the Aegean extensional province of western Turkey (e.g., McKenzie, 1972). The fault surface is exposed by erosion and mining activity southeast of the town of Kuşadasi, Aydin District, in western Turkey (37°49′29″N; 27°15′57″E). The surface model is an ∼8 m × 5 m area of an outcrop in Tertiary limestone described in detail by Hancock and Barka (1987).

In the model area, the fault surface is characterized by fault wear surface features, such as long-wavelength, slip-parallel corrugations plunging steeply toward the left (southwest) (Fig. 5B), and comb fracture traces (Hancock and Barka, 1987) oriented at a high angle to the corrugations Figs. 5A, 5B). The maximum roughness at a compass length range between 2 cm and 0.35 mm was measured parallel to slip, and perpendicular to the comb fractures, and the anisotropy is moderate with *a* = 0.69 (Fig. 5C). Smaller scale features produce a much less coherent roughness plot with a much smaller anisotropy of *a* = 0.84 (Fig. 5D), possibly including the scratch marks of excavation machinery (jagged white lines in Fig. 5A). This result is different from fault surface analyses carried out by Power and Tullis (1991) and Sagy et al. (2007), who found that the surfaces were smoothest parallel to the direction of slip. In our image the comb fractures, which are oriented at a high angle to fault slip, dominate the roughness, while the wavelength of the fault surface corrugation, in the meter range, is too large to be detected.

#### Heavitree Quartzite at Ormiston Gorge

The Neoproterozoic Heavitree quartzite is characterized by abundant joints and veins. The joints display a very regular and near “Cartesian” geometry (Hobbs, 1993), which may be related to folding and thrusting (e.g., Flöttmann and Hand, 1999). The surface model we analyzed is an ∼2 m × 1.4 m area of a gorge wall, which has been polished by the episodically flowing Finke River in Ormiston Gorge, Northern Territory, Australia (23°37′38″S; 132°43′38″E).

The rock face morphology in the model is dominated by a set of fractures dipping steeply toward the left (∼south), and a subhorizontal set. The roughness plot at compass lengths between 2 cm and 0.35 mm detects a maximum at ∼129°, which may be a composite value (Fig. 6C); anisotropy is moderate at 0.72. At a compass length range between 0.35 and 0.06 mm the roughness is dominated by a maximum at ∼24°, which corresponds well with the perpendicular orientation of a small vein set (Fig. 6A). Roughness anisotropy at this scale is higher with *a* = 0.62 (Fig. 6D). An ∼20 cm surface offset along a steep fracture in the middle of the image (Figs. 6A, 6B) is not detected because its amplitude is too large compared to the compass lengths.

#### Taunus Quartzite at Trechtingshausen Quarry

The Lower Devonian Taunus quartzite is mined for aggregate in a large quarry in Trechtingshausen, Rhineland-Palatinate, Germany (50°01′00″N; 7°49′45″E). The model is ∼10 m × 6.5 m, and represents a fresh vertical surface shaped by blasting (Fig. 7A). The rock fabric in the Trechtingshausen pit is dominated by a number of discontinuity plane sets, including shallowly dipping bedding planes, a horizontal cleavage set, and a near vertical fracture cleavage. The planar deformation fabrics formed during movement along the Taunuskammüber-schiebung, a major Variscan thrust zone in the Rhenohercynian slate belt (e.g., Oncken, 1988), as well as during the Quaternary uplift of the Rhenish Massif (e.g., Meyer and Stets, 1996). At compass lengths between 64 cm and 11 cm (Fig. 7C) these fabric elements generate roughness maxima oriented ∼135° (perpendicular to bedding), and ∼35° (perpendicular to a fracture set). There is a general dominance of horizontal roughness, probably related to the near vertical fracture set. Anisotropy is moderate at 0.72. At a compass length range between 11 and 2 cm roughness is dominantly subparallel to bedding, which dips ∼30° to the left in Figure 7A. Anisotropy is much lower at *a* = 0.83 (Fig. 7D).

## DISCUSSION

We have developed a technique to detect anisotropic features on rock faces based on fractal analysis. Our preliminary results from a variety of natural examples suggest that there can be a significant change of roughness value, of anisotropy, and of anisotropy direction, when changing the scale of observation, which confirms earlier research proposing a self-affine rather than a self-similar fractal geometry for rock surface characteristics (i.e., Kulatilake and Um, 1999, and references therein). The method we have developed is an important step toward developing a toolkit to automatically detect and interpret mechanical and chemical rock surface characteristics from digitally acquired data sets. Our examples of roughness measurements are based on data acquired by digital photogrammetry, but the analysis can be applied to point cloud representation of any surface. Future improvements to the method include techniques to detect how roughness varies spatially across subdomains of the rock face, as well as using Fourier series and wavelet analysis techniques to calculate industry standards, such as the joint roughness coefficient (JRC). The application of the method can be extended to systematic study of rock surface features across scale, and to surface analysis of geometrical objects with closed surfaces in three dimensions.

When applying this approach to 3D point clouds of surfaces, careful consideration must be given to how the scale of roughness measurement relates to the spatial resolution of the input data. In our experiment, we used an average distance between adjacent points within the point cloud as a guide to determine our minimum compass length. As the profiles are sampled at various locations, the profile lengths vary. This needs to be taken into account when selecting the maximum compass length. In our analysis we have chosen an arbitrary value of ∼1/20 of the diagonal length across the 3D model as a guide to determine the maximum compass lengths. The error within the fractal dimension caused by a selection of suboptimal compass length ranges can be reduced by using fractal interpolation methods (Xie et al., 2001; Gemperline and Siller, 2002). To ensure correct results the precision of the point cloud should be taken into account. The precision of the point cloud will ultimately depend on the photogrammetry workflow, particularly factors like the resolution of the camera, and accuracy of the surveyed points used to reference the model.

There are many fractal analysis methods that are previously applied for profile analysis. Schmittbuhl and Vilotte (1995), for example, discussed various methods applied to synthetically generated self-affine profiles. Our experiments indicate that the error inherent to the compass method does not affect the detection of the directions of minimum and maximum roughness.

## CONCLUSIONS

We have developed a method of measuring the variation of surface roughness across different orientations of a surface. This method detects the directions of minimum and maximum roughness, and gives an indication on which scale of observation the anisotropy is most pronounced. This measurement technique has a potential for widespread application in rock engineering and structural geology.

We thank Michael Drews and Hagen Deckert for supplying the Taunus quartzite dataset; Cameron D. Walsh for preprocessing of the Cullercoats images; Bob Holdsworth, Jonny Imber, Steven Smith, Adam Pugh, Ruth Wightman, Juilin Guo, Robert Wilson, Ken McCaffrey, and Richard Jones for organizing the Penrose conference field trip to Cullercoats; George Poropat and Wayne Robertson for help and support with Sirovision; Talip Güngör for assistance in Turkey; Paul Bourke for help with fractal dimension calculations; and iVEC ‘The hub of advanced computing in Western Australia’ for funding Baker's internship. We also thank Steven Smith and an anonymous reviewer for constructive and very helpful suggestions.