Box-counting methods to directly estimate the fractal dimension of a rock surface

Applied Surface Science 314 (2014) 610–621

Contents lists available at ScienceDirect

Applied Surface Science journal homepage: www.elsevier.com/locate/apsusc

Box-counting methods to directly estimate the fractal dimension of a rock surface T. Ai a , R. Zhang b,∗ , H.W. Zhou c , J.L. Pei b a

School of Architecture and Environment, Sichuan University, Chengdu 610065, China State Key Laboratory of Hydraulics and Mountain River Engineering, College of Hydraulic and Hydropower Engineering, Sichuan University, Chengdu 610065, China c School of Mechanics, Architecture and Civil Engineering, China University of Mining and Technology, Beijing 100083, China b

a r t i c l e

i n f o

Article history: Received 9 December 2013 Received in revised form 10 June 2014 Accepted 13 June 2014 Available online 11 July 2014 Keywords: Rock surface Box-counting dimension Cubic covering Takagi surface

a b s t r a c t Surfaces of rocks are usually not perfectly “smooth”, and two box-counting methods, i.e. the conventional cubic covering method (CCM) and improved cubic covering method (ICCM), can directly describe the irregularities of a rock fracture surface without any approximate calculations. Our investigation showed that if the scale ı of covering cubes is greater than the sampling interval S0 , the CCM and ICCM cannot completely cover the object rough surface. Considering this, we presented two new cubic covering methods, namely the differential cubic covering method (DCCM) and relative differential cubic covering method (RDCCM) to directly evaluate the fractal dimension of a rough surface according to the definition of box-counting dimension. Experimentally, a 3D laser profilometer was used to measure the topography of a natural surface of sandstone. With the CCM, ICCM, DCCM and RDCCM, direct estimations of the fractal dimension of the rock surface were performed. It was found the DCCM and RDCCM usually need more cubes to cover the whole fracture surface than the CCM and ICCM do. However, the estimated fractal dimensions by the four methods were quite close. Hence, three Takagi surfaces with known fractal dimensions of 2.10, 2.50 and 2.90 were adopted to further examine the four boxcounting algorithms. Results showed that for a low fractal dimension Takagi surface, the DCCM and RDCCM gave accurate results within the ranges determined by small covering scales, whereas the CCM and ICCM always overestimate the fractal dimension for all the potential scale ranges investigated in current work; for high fractal dimension surfaces, the CCM and ICCM provided very good results within the ranges determined by small covering scales, and oppositely, the DCCM and RDCCM cannot provide a good estimation of the fractal dimension within such scale ranges but can determine approximate results at large scales. © 2014 Elsevier B.V. All rights reserved.

1. Introduction The classic methods of solid mechanics assume the homeomorphism of deformation [1]. However, voids, cracks, faults, fractures and joints always exist in rocks as a result of various geological processes [2]. The mechanical properties of fractured rocks principally depend on the state of existing discontinuities [3], and moreover, the behaviors of these discontinuities are strongly affected by their surface characteristics [4]. Therefore, the quantitative description of the topography of a fracture surface is very important.

∗ Corresponding author. Tel.: +86 02885465866. E-mail addresses: [email protected] (T. Ai), [email protected] (R. Zhang), [email protected] (H.W. Zhou), [email protected] (J.L. Pei). http://dx.doi.org/10.1016/j.apsusc.2014.06.152 0169-4332/© 2014 Elsevier B.V. All rights reserved.

Material surfaces, especially the surfaces of rock-like materials, are usually not perfectly “smooth” and the irregularities in the form of valleys and convexities always exist [5]. Based on the fractal geometry [6], a theory to characterize the degree of irregularities, many researchers quantitatively described the morphology of rock fracture surfaces, and it has been confirmed that fracture surfaces in rocks exhibited a statistical fractal behavior in a certain scale range [7]. The fractal dimension D can be used to measure the irregularities and the degree of complexity of surface shape [5] and the intercept, in a log–log way, is an indicator of asperities [8]. In the early stage, in order to simplify the problem and avoid data acquisition difficulties, a linear sectional profile of a surface was widely characterized to grasp the roughness of a 3D surface [9]. The indirect measurement methods (slit island analysis, the divider and the self-affine variogram) were often employed to measure such sectional fracture profiles [10]. As a

T. Ai et al. / Applied Surface Science 314 (2014) 610–621

result, fractal dimensions measured by these methods fall in the range 1 < D < 2 [10]. However, a one-dimensional analysis provides an incomplete and even biased characterization of a fracture surface [11]. Thus, the fractal dimension of a rough surface obtained by adding 1.0 to the fractal dimension obtained from a single sectional profile of that surface was widely adopted [11]. Such an approximation might be very close to the real fractal dimension of a fracture surface, but essentially, it is unacceptable [11]. In order to perform a direct measurement of the fractal dimensions of fracture surfaces, the measurement technique should be taken into account first [10]. As summarized in [12,13], the scanning electron microscopy (SEM) [14,15], atomic force microscopy (AFM) [16–19], mechanical stylus profilometry (MSP) [20] and non-contact laser profilometry (LPM) [21,22] are commonly used to conduct a topographical measure of a surface. Among these methods, the profilometric methods can provide quantitative topographical information of a rough surface [12,13]. Especially the non-contact optical instruments, such as a 3D laser profilometer, make it more convenient to obtain the morphology data of a rock fracture surface in a form of (x, y, z) coordinates. In the theoretical aspect, the triangular prism surface area method (TPSAM) [23], proposed in 1986, made it possible to conduct a direct measurement of the fractal dimensions of a fracture surface. In this method, a three-dimensional geometric equivalent of the “walking dividers” method in two dimensions was adopted [23]. Elevation values at the corners of cubes were employed to interpolate a center value. However, the elevation at the center of each grid cell is determined by linear interpolation of the four heights of the adjacent points. Thus, it is almost impossible to exactly calculate the true area of the fracture surface within the grid cell [11]. Considering this, the projective covering method (PCM) [10], which provides a similar fractal dimension result, was proposed. In this method, the real area surrounded by four points on the fracture surface is approximated by two triangles and every point for calculation of the approximate area can be assured to be on the fracture surface. Recently, Kwa´sny [24] modified the PCM by introducing a more precise area calculation method. Nevertheless, the area of the fracture surface is also approximate. To accurately estimate the fractal dimension of a fracture, a new method, namely cubic covering method (CCM) [11] was presented according to the principle of the covering method. The most important consideration is that every point (which is essentially the sampling point) used is exactly located on the rough surface, thus the calculated fractal dimension is absolutely accurate. Based on the CCM, an improved cubic covering method (ICCM) was proposed [25]. The cubic covering procedure was implemented from a universal basis plane for each grid with a measurement scale ı for the ICCM. Generally, for a box-counting method, it should be able to completely cover a fractal set [26,27]. However, it is almost impossible to ensure that the maximum (minimum) one among the height values of the four intersection points is exactly the maximum (minimum) height of the irregular surface area within the scale ı when ı is larger than the sampling interval of the laser scanning due to the complexity of surface shape. Therefore, both the CCM and ICCM, strictly speaking, cannot totally cover a fractal set. Considering this, we presented two new cubic covering methods, namely the differential cubic covering method (DCCM) and relative differential cubic covering method (RDCCM). Both the DCCM and RDCCM can totally cover the fracture surface. Later, a laser profilometer was then used to obtain the data set of a fracture surface of sandstone. With the data set, we directly determined the fractal dimension of the fracture surface of sandstone by the CCM, ICCM, DCCM and RDCCM, respectively, and quite similar results were obtained. Finally, series of rough surfaces which have known fractal dimensions were generated based on

611

the Takagi function to validate the applicability of the four covering methods. 2. Box-counting algorithms 2.1. Box-counting dimension Due to the ease of mathematical calculation and empirical estimation, box-counting dimension is one of the most widely used dimension [27]. To find the box-counting dimension of F, a nonempty bounded subset of Rn , one may draw a mesh of squares or boxes of side ı and count the number N(ı) that overlap the set [27]. Through changing the scale ı, different values of N(ı) can be obtained. The total number N(ı) of cubes depends on the used measurement scale, ı. If the fracture surface exhibits the fractal behavior, the relation between N(ı) and ı is given by N(ı)∼ı−D ,

(1)

where D is the fractal dimension of the object fracture surface. Log transformation of this simple power law yields a straight line with slope −D. Therefore, the box-counting dimension can be estimated by the gradient of the graph of ln N(ı) against ln ı given by DB = −lim

ı→0

ln Nı (F) . ln ı

(2)

The number N(ı) of ı-mesh cubes that intersect a set is an indicator how spread out or irregular the set is when examined at a scale of ı [27]. The box-counting dimension shows how rapidly the irregularities develop as ı tends to 0 [27]. 2.2. The conventional cubic covering method and its improved version The conventional cubic covering method [11] provides a very simple way to directly calculate the fractal dimension of a rough surface. It is assumed that there exists a regular cube grid on the plane XOY (see Fig. 1a), and that in each grid cell with scale ı, four intersection points correspond to four heights of a fracture surface: h1 (i, j), h2 (i, j + 1), h3 (i + 1, j), and h4 (i + 1, j + 1) (where 1 ≤ i, j ≤ n − 1, n is the total number of sampling points on each individual profile on a fracture surface). If the cubes with the measurement scale ı are used to cover the irregular surface area within the scale ı, the maximum difference among h1 (i, j), h2 (i, j + 1), h3 (i + 1, j), and h4 (i + 1, j + 1) will determine the number Ni,j of the required cubes: Ni,j = INT{ı−1 [max(h1 (i, j), h2 (i, j + 1), h3 (i + 1, j), h4 (i + 1, j + 1)) − min(h1 (i, j), h2 (i, j + 1), h3 (i + 1, j), h4 (i + 1, j + 1))] + 1},

(3)

where INT rounds the element to the nearest integer towards positive infinity. Then, the total number of cubes required for covering the whole fracture surface is N(ı) =

n−1 

Ni,j .

(4)

i,j=1

Apparently, the covering process in each grid cell is always implemented from the lowest point among the four intersection points. Zhang et al. [25] pointed out that it is difficult to well describe the complexity of a rough surface with such a covering procedure, and therefore, they proposed an improved cubic covering method (ICCM, the schematic view was shown in Fig. 1b). The ICCM gives the number Ni,j of cubes need to cover the irregular surface area within the scale ı by

612

T. Ai et al. / Applied Surface Science 314 (2014) 610–621

Ni,j = INT[ı−1 max(h1 (i, j), h2 (i, j + 1), h3 (i + 1, j), h4 (i + 1, j + 1))] − INT[ı−1 min(h1 (i, j), h2 (i, j + 1), h3 (i + 1, j), h4 (i + 1, j + 1))] + 1.

(5)

In order to better show the covering processes of the CCM and ICCM, a hypothetic linear profile (denoted by the red solid line in Fig. 2) was covered with the reduced forms of the two methods. The sample profile was determined by “hypothetically scanning” with a sampling interval S0 and nine sampling points, respectively denoted by the letters A, B, . . ., I, were obtained. The maximum height among the sampling points was 8S0 . The scale of covering cubes was set to S0 , 2S0 and 4S0 for both covering methods. For a 2D case, Eqs. (3) and (5) can be respectively simplified as Ni = INT{ı−1 [max(h1 (i), h2 (i + 1)) − min(h1 (i), h2 (i + 1))] + 1}, (6) and Ni = INT[ı−1 max(h1 (i), h2 (i + 1))] −INT[ı−1 min(h1 (i), h2 (i + 1))] + 1.

(7)

As shown in Fig. 2, if ı = S0 , the CCM covered the segment AB from point A and four cubes were required according to Eq. (6); the ICCM covered AB from the basis plane (x-axis in this example) and four cubes were required according to Eq. (7). For the entire profile, 28 and 31 cubes were respectively required for the CCM and ICCM. If ı = 2S0 , the CCM implemented the covering procedure from point E to cover the segments CD and DE and according to Eq. (6), two cubes were required; the ICCM used three cubes to cover the segments CD and DE from the basis plane. Eight and eleven cubes were respectively required for the CCM and ICCM to cover the sample profile. When ı = 4S0 , the CCM started to cover the segments EF, FG, GH and HI from point I, and according to Eq. (6), only one cube was required. According to Eq. (7), the ICCM also employed one cube to cover these four segments. For both of the CCM and ICCM, two cubes were adopted to cover the entire profile. 2.3. The differential and relative differential cubic covering methods

Fig. 1. Schematic views of the covering methods. (a) CCM (according to [11]), (b) ICCM, (c) DCCM, (d) RDCCM.

From the above descriptions, it is clear that if ı = 2S0 , the heights of the points A, C, E, G and I determined the number Ni ; if ı = 4S0 , only three points, A, E and I, were used to determine Ni . In other words, if the scale ı of covering cubes is greater than the sampling interval S0 , the total number of the cubes needed to cover the entire profile or fracture surface depends on some specific sampling points. According to definition of box-counting dimension, the cubes that overlap or intersect the set F should be counted [27]. Under this consideration, an important issue is to keep all the measured information of the surface. As shown in Fig. 1c, four intersection points and [(ı/S0 + 1)2 − 4] intermediate points correspond to (ı/S0 + 1)2 heights of the fracture surface area in the field of the (i, j)th cell unit: hk,q (i, j) (where k and q are the sub-numbers of a sampling point on the fractufe surface area; 1 ≤ k, q ≤ ı/S0 + 1; i and j are the numbers of a grid cell with the scale ı on each individual profile direction of a fracture surface; 1 ≤ i ≤ (n – 1) × S0 /ı, 1 ≤ j ≤ (m – 1) × S0 /ı, m and n are the total numbers of sampling points on each individual profile direction). If the cubes with scale ı are used to cover the irregular surface area from the reference

T. Ai et al. / Applied Surface Science 314 (2014) 610–621

613

Fig. 2. Differences of covering process of the CCM, ICCM, DCCM and RDCCM. (For interpretation of the references to color in the text citation to this figure, the reader is referred to the web version of this article.)

614

T. Ai et al. / Applied Surface Science 314 (2014) 610–621

plane, the maximum and minimum heights of the intersection and intermediate points will determine the number Ni,j of cubes needed to cover the irregular surface area in the field of the (i, j)th cell unit, i.e. Ni,j = INT

1 ı

− INT





max h1,1 (i, j), ..., hk,q (i, j), ..., hı/S0 +1,ı/S0 +1 (i, j)

1





min h1,1 (i, j), ..., hk,q (i, j), ..., hı/S0 +1,ı/S0 +1 (i, j)

ı

+ 1.

(8)

Then, the total number of cubes needed to cover the whole fracture surface is

⎡



(n−1)·S0 /ı

N(ı) =



(m−1)·S0 /ı

⎣

i=1

⎤

Ni,j ⎦ .

(9)

j=1

If the fracture surface exhibits fractal behavior, the fractal dimension of the rough surface is also given by Eq. (1). The related algorithm represented by Eqs. (8) and (9) is named as the differential cubic covering method (DCCM). In addition, one can also cover the fracture surface area in the field of the (i, j)th cell unit from the lowest point of the surface area (see Fig. 1d). In such a case, the relative difference between the maximum and minimum values of heights of the surface area will determine Ni ,j , i.e. Ni,j = INT

1





max h1,1 (i, j), ..., hk,q (i, j), ..., hı/S0 +1,ı/S0 +1 (i, j)

ı





−min h1,1 (i, j), ..., hk,q (i, j), ..., hı/S0 +1,ı/S0 +1 (i, j)



+1 . (10)

The total number of cubes need to cover the whole fracture surface is also determined by Eq. (9), and the fractal dimension of the rough surface is given by Eq. (1). This box-counting algorithm was called as the relative differential cubic covering method (RDCCM). In order to better understand the DCCM and RDCCM, we also take the linear profile in Fig. 2 as a sample profile. For a 2D fractal analysis, Eqs. (8) and (10) can be respectively reduced as Ni = INT

1 ı

− INT and Ni = INT





max h1 (i), ..., hk (i), ..., hı/S0 +1 (i)

1 ı

1 ı









min h1 (i), ..., hk (i), ..., hı/S0 +1 (i)

+ 1,

(11)

max h1 (i), ..., hk (i), ..., hı/S0 +1 (i)





− min h1 (i), ..., hk (i), ..., hı/S0 +1 (i)



+1 .

(12)

And Eq. (9) can be simplified as



(n−1)·S0 /ı

N(ı) =

Ni .

(13)

i=1

When ı = S0 , no intermediate point exists and Eqs. (11) and (12) can be written as Ni = INT

1 ı



max(h1 (i), h2 (i)) − INT

1 ı



min(h1 (i), h2 (i)) + 1, (14)

and Ni = INT

1 ı



[max(h1 (i), h2 (i)) − min(h1 (i), h2 (i))] + 1

,

(15)

respectively. In this case, Eqs. (7) and (14) are essentially the same, and Eqs. (6) and (15) are exactly the same. Therefore, both the DCCM and ICCM have the same covering results for ı = S0 (see Fig. 2d and g). Similarly, as shown in Fig. 2a and j both the CCM and RDCCM gave the same result. When ı = 2S0 , there are four intermediate points B, D, F and H. For the segments EF and FG, the CCM started the covering procedure from point G and two cubes were required according to Eq. (6); the ICCM started to cover the same segments from x-axis and three cubes were needed according to Eq. (7). Actually, the point F was the highest one among the points E, F and G. According to the presented DCCM and RDCCM, the required cubes could be determined by the highest point (F) and the lowest one (G). According to Eqs. (11) and (12), four cubes were actually required for either of the DCCM and RDCCM. Additionally, from Fig. 2a and j, it can be seen that the DCCM and RDCCM needed three and two cubes respectively to cover the segments CD and DE. When ı = 4S0 , the points B, C, D, F, G and H were the intermediate points. In the CCM and ICCM, the points A, E and I were used to calculate the number of the cubes needed to cover the whole profile, and the number was two for the two methods. For the segments AB, BC, CD and DE, the highest point (the intermediate point C) and the lowest point (the intersection point A) were adopted to determine the number of the covering cubes by the DCCM and RDCCM. Similarly, the highest intermediate point F and the lowest intermediate point G were employed to calculate the number of the cubes needed to cover the segments EF, FG, GH and HI, and the calculated number was two. Therefore, to cover the whole profile, both the DCCM and RDCCM needed four cubes. Obviously, the DCCM and RDCCM kept all the height information of the measured rough profile or fracture surface. 3. Fractal analysis of a real rock surface 3.1. Surface topography measurement In this paper, the testing material, sandstone, was acquired from a stone quarry in Fangshan County, Beijing. Its dry density is about 2.41 g/cm3 . The sandstone core was drilled with a diamond drill bit with a diameter of 50 mm. With a rotary diamond saw, a circular disk was then cut from the core and the thickness of the disk was approximately 25 mm. A fresh and unfilled fracture surface was produced in the sample with the Brazilian method. Topographical measurement of the fracture surface was carried out with a 3D laser profilometer which is a noncontact device and enables one to obtain data files of x, y and z coordinates. Its laser probe provides an accuracy of ±0.1 mm and an elevation range of 300 mm. It can move automatically over the sample according to the programmed path to measure the topography of the fracture surface. The obtained data set of a fracture surface topography consists of both coordinates (x, y) and corresponding heights (z) of the object surface. The fracture surface was scanned with a sampling interval of 0.1 mm in this study. Then, an isometric view of an internal scanning field (24 mm × 48 mm) was reconstructed (as shown in Fig. 3). 3.2. Fractal analysis by the CCM, ICCM, DCCM and RDCCM The data set of the topography of the sandstone surface (see in Fig. 3) was employed to make direct determinations of the fractal dimension with the CCM, ICCM, DCCM and RDCCM. For each method, the relations between N(ı) and ı are given in Table 1. It

T. Ai et al. / Applied Surface Science 314 (2014) 610–621

Fig. 3. Surface topography of a rock fracture measured by the laser profilometer (scanning interval = 0.1 mm; total number of measured points = 240 × 480).

Table 1 Relations between N(ı) and ı estimated by the CCM, ICCM, DCCM and RDCCM. ı

N(ı) CCM 0.1 0.2 0.3 0.4 0.5 0.6 0.8 1.0 1.2 1.5 1.6 2.0 2.4 3.0 4.0 4.8 6.0 8.0 12.0 24.0

195,394 41,925 16,862 9100 5555 3718 2047 1254 838 529 459 289 201 128 72 50 32 18 8 2

ICCM 251,913 55,687 22,716 12,294 7528 5006 2784 1658 1135 720 641 393 249 152 100 63 32 18 8 2

DCCM 251,913 59,765 25,223 13,974 8645 5766 3245 2000 1337 845 738 446 301 191 115 67 32 18 8 2

RDCCM 195,394 44,841 18,549 10,039 6211 4191 2287 1456 999 598 517 300 202 128 72 50 32 18 8 2

can be seen that within the scale range of 0.1–4.8 mm, the numbers of required cubes for each method have the following order: NDCCM (ı) > NICCM (ı) > NRDCCM (ı) > NCCM (ı), where the footnotes are the names of each box-counting method, while their results are the same within the scale range of 6.0–24.0 mm. When the covering scale ı equals the sampling interval S0 , 195,394 cubes were needed for both the CCM and ICCM, and 251,913 cubes for both the DCCM and RDCCM. Fig. 4 shows the log–log plots of N(ı) and ı estimated by each method. For the CCM, the data points fall into two straight segments. The correlation coefficients of linear regression analysis are respectively 0.9998 on the scale ranging from 0.1 to 2.4 mm, and 1.0 on the scale ranging from 3.0 to 24.0 mm. As shown in Fig. 4a, within the scale range of 0.1–2.4 mm, the fractal dimension of the fracture surface is 2.1655, and on the scale ranging from 3.0 to 24.0 mm, the fractal dimension is 2.0. Obviously, within the lower measurement resolution range from 3.0 to 24.0 mm, the fracture surface does not exhibit fractal behavior at all. When the ICCM was adopted to analysis the rock surface, the points denoted by the black solid symbols

615

at the measurement scales of 4.0 mm and 4.8 mm deviate from the regression lines determined by the points within the scale ranges of 0.1–3.0 mm and 6.0–24.0 mm (as shown in Fig. 4b) and the data points fall into three segments. On the scale ranging from 0.1 to 3.0 mm, the correlation coefficient of linear regression analysis and the fractal dimension are 0.9999 and 2.1662, respectively; on the scale ranging from 6.0 to 24.0 mm, the correlation coefficient and the fractal dimension are 1.0 and 2.0, respectively. As shown in Fig. 4c, for the DCCM, one “singular” point was observed to deviate from the regression lines respectively determined by points within the scale ranges of 0.1–4.0 mm and 6.0–24.0 mm. The fractal dimension of the fracture surface is 2.1084 on the scale ranging from 0.1 to 4.0 mm; within the scale range of 6.0–24.0 mm, the fractal dimension is 2.0. The data points fall into two straight segments if the RDCCM was employed, see in Fig. 4d. The correlation coefficients of linear regression analysis are respectively 0.9998 on the scale ranging from 0.1 to 2.4 mm, and 1.0 on the scale ranging from 3.0 to 24.0 mm, and the fractal dimensions for each scale range are 2.1524 and 2.0, respectively. Additionally, from Fig. 4, it can be seen that the fracture surface exhibits fractal behavior within the measurement scale range of 0.1–2.4 mm for the CCM and RDCCM, range of 0.1–3.0 mm for the ICCM and range of 0.1–4.0 mm for the DCCM. In Fig. 4, it is obvious that the fractal dimensions at the smaller scales are higher than those at the larger scales for all the methods. Hence, in order to obtain an accuracy estimation of the roughness of a fracture surface, smaller covering scales are required, and a higher topographical measurement resolution is also required for reflecting the influence of the structures of the fracture surface on fractal characteristics [11]. 4. Uncertainty in fractal dimension estimation by the CCM, ICCM, DCCM and RDCCM Compared with the TPSAM, PCM and the modified PCM, the CCM, ICCM, DCCM and RDCCM can avoid the approximate estimation of the real area surrounded by four points on the fracture surface and assure that every covering step is accurate. However, the results estimated by the CCM, ICCM, DCCM and RDCCM in the previous section differed only by about 2.71%. Therefore, it may fall within the uncertainty of the applied methods and a discussion on this issue is necessary. 4.1. Rough surface modeling with Takagi function To validate the box-counting methods, series of Takagi surfaces with known fractal dimensions were generated. Takagi functions are continuous nowhere differentiable and can be seen as a superposition of “pyramids” of various frequencies attenuated by a scaling factor [28]: z(x, y) =

∞ 

bn (2n−1 x, 2n−1 y),

(16)

n=1

where 1/2 < b < 1 and (x, y) ∈ I2 . The function , called the generating kernel, is a separable function that consists of the product of two triangular functions [28], i.e. (x, y) = |2x − fix(2x)| · |2y − fix(2y)|,

(17)

where fix(x) rounds x to the nearest integer towards zero. The fractal dimension of a Takagi surface is between 2 and 3, and is given by DTakagi =

log(8b) . log 2

(18)

In a real surface modeling case, the value of n in Eq. (16) cannot be infinite. The more accurate results of z(x, y), the larger the n

616

T. Ai et al. / Applied Surface Science 314 (2014) 610–621

6.0

6.0 5.0

log N = -2.1655 log(δ) + 3.1039 R² = 0.9998

4.0 log [N(δ)]

log [N(δ)]

4.0 3.0

3.0 2.0

2.0

1.0

1.0

log N = -2 log(δ) + 3.0615 R² = 1 -1.2 -0.9 -0.6 -0.3

0.0 0.3 0.6 log (δ) / mm

0.9

1.2

-1.2 -0.9 -0.6 -0.3

1.5

(a) Covering method: CCM

0.9

1.2

1.5

1.2

1.5

6.0 5.0

log N = -2.1084 log(δ) + 3.2984 R² = 0.9999

log N = -2.1524log(δ) + 3.1465 R² = 0.9998

4.0

log [N(δ)]

4.0 log [N(δ)]

0.0 0.3 0.6 log (δ) / mm

(b) Covering method: ICCM

6.0

3.0

3.0

2.0

2.0 1.0

log N = -2 log(δ) + 3.0615 R² = 1

0.0

0.0

5.0

log N = -2.1662 log(δ) + 3.2294 R² = 0.9999

5.0

1.0

log N = -2 log(δ) + 3.0615 R² = 1

log N = -2log(δ) + 3.0615 R² = 1

0.0

0.0 -1.2 -0.9 -0.6 -0.3

0.0 0.3 0.6 log (δ) / mm

0.9

1.2

1.5

(c) Covering method: DCCM

-1.2 -0.9 -0.6 -0.3

0.0 0.3 0.6 log (δ) / mm

0.9

(d) Covering method: RDCCM

Fig. 4. Log–log plots of N(ı) and ı estimated by the CCM, ICCM, DCCM and RDCCM.

used. However, it would make little differences in the results of z(x, y) after n exceeds a certain value. In this work, after checking the results with n = 20 and n = 50, no obvious variations were found. Therefore, only the first 50 elements of the series in Eq. (16) were used. The fractal surfaces with fractal dimensions 2.10, 2.50 and 2.90 generated by using nmax = 50 were presented in Fig. 5, and in all the surface generations the net of 64 × 64 squares was used.

4.2. Applicability of the applied box-counting algorithms The Takagi surfaces showed in Fig. 5 were employed to validate the applicability of the CCM, ICCM, DCCM and RDCCM. To guarantee the calculation accuracy, in all surface generation cases the net of 215 × 215 squares was used. That is saying a very small virtual scanning interval S0 of 1/(215 ) was adopted to virtually scan

Fig. 5. Takagi surfaces with different fractal dimensions.

Table 2 Fractal dimensions estimated by each box-counting method for full scale ranges, whereas the theoretical value is exactly 2.10.

T. Ai et al. / Applied Surface Science 314 (2014) 610–621 617

618

Table 3 Fractal dimensions estimated by each box-counting method for full scale ranges, whereas the theoretical value is exactly 2.50.

T. Ai et al. / Applied Surface Science 314 (2014) 610–621

Table 4 Fractal dimensions estimated by each box-counting method for full scale ranges, whereas the theoretical value is exactly 2.90.

T. Ai et al. / Applied Surface Science 314 (2014) 610–621 619

620

T. Ai et al. / Applied Surface Science 314 (2014) 610–621

all the rough surfaces. Therefore, points with a total number of 1,073,807,361 were obtained for each Takagi surface. The covering scales ı (lengths of the cube sides) were set as: ıi = 2(i−1) × S0 (i = 1, . . ., 14). The results of a real rock surface showed in Fig. 4 indicated that there is no universal fractal dimension on all scales, which is consistent with the results in [11,22,25,29]. Obviously, the scale range significantly influences the estimated result of fractal dimension. Therefore, we used all the potential scale ranges, e.g. range from ı1 to ı2 , range from ı4 to ı13 and etc., to better examine the scale range effects and the applicability of each box-counting method. Fractal dimension results estimated by the CCM, ICCM, DCCM and RDCCM for the three Takagi surfaces within all potential scale ranges are presented in Tables 2–4. In these tables, the fractal dimension in the grid overlapped by the row corresponded to ıi (i = 1, . . ., 13) and the column corresponded to ıj (j = 2, . . ., 14) was determined by the scale range from ıi to ıj . For example, in Table 2c, the result of 2.35 was obtained within the range ı11 to ı14 by taking log transformation of Eq. (1) to yield a straight line with slope −2.35. In Table 2, the fractal dimension of the adopted Takagi surface is exactly 2.10. However, it was found that the CCM and ICCM cannot provide a good estimation of the fractal dimension for all the studied scale ranges in this work, whereas the DCCM and RDCCM can give very accurate results within the ranges determined by small covering scales (e.g. within the range ı1 to ı6 , both the DCCM and RDCCM gave a result of 2.10). For Table 3, the fractal dimension of the adopted Takagi surface is 2.50. It was clear that the DCCM and RDCCM can well estimate the fractal dimension of the rough surface at large scales and the CCM and ICCM can give very good results within the ranges determined by small covering scales (e.g. within the range ı1 to ı6 , both the CCM and ICCM gave a very close result of 2.51; within the range ı5 to ı14 , the DCCM provided a result of 2.51; and within the range ı4 to ı14 , the RDCCM gave a value of 2.50). In Table 4, the fractal dimension of the Takagi surface is 2.90, and the CCM and ICCM also gave very good results at small covering scales, and oppositely, the DCCM and RDCCM cannot provide a good estimation of the fractal dimension at such scales. Both the DCCM and RDCCM can determine approximate results compared to the theoretical fractal dimension within the ranges determined by larger scales. Additionally, from Tables 2–4, it can be seen that the DCCM and RDCCM gave the fractal dimension values in a range of 2.0–3.0 for the three surfaces, whereas the CCM and ICCM might give fractal dimension values larger than 3.0 within some scale ranges for surfaces with high fractal dimensions. From the above descriptions, it can be concluded that the uncertainty in fractal dimension estimation by the CCM, ICCM, DCCM and RDCCM exists. Within a range determined by small coverig scales, the DCCM and RDCCM gave accurate results for low fractal dimension Takagi surfaces, whereas the CCM and ICCM provided very good results for high fractal dimension Takagi surfaces.

5. Discussion and conclusion In many fields of rock mechanics, both of the characterization and re-modeling of a rough surface are very important. Under this consideration, the data set consists of the x, y and z coordinates of a rough surface which can enable one to directly characterize the morphology and conveniently re-model the natural rock fracture is disabled. Using such a coordinate data set, the CCM, ICCM, DCCM and RDCCM can directly and easily estimate the fractal dimension. According to the fractal analysis results of a sandstone fracture surface topography, it was found within the covering scale ranges determined by large scales, e.g. 3.0–24.0 mm for the CCM and RDCCM and 6.0–24.0 mm for the ICCM and DCCM, the fracture surface does not exhibit fractal behavior. Hence, in order to

obtain an accuracy estimation of the roughness of a fracture surface, smaller covering scales are required. However, it seems completely opposite if a Takagi surface was employed to conduct a fractal analysis. From Tables 2–4, it can be seen a higher fractal dimension will be obtained if a scale range combined by large measurement scales is used. One possible explanation is that the natural fractals are more complicated than a theoretical or mathematical fractal. However, such an interesting subject needs more investigations. In this paper, the CCM and ICCM were found difficult to totally cover a fractal set. Considering this, two new cubic covering methods, i.e. the DCCM and RDCCM were presented in this paper and the two methods can completely cover the fracture surface. Then, the four box-counting methods were used to analyze a rough surface of sandstone and three Takagi surfaces with known fractal dimensions. Results showed that: (1) for a natural rock surface of sandstone, there is no universal fractal dimension at all scales and the estimated fractal dimensions by these four methods were quite close; the natural sandstone surface exhibits fractal behavior within the measurement scale range of 0.1–2.4 mm for the CCM and RDCCM, range of 0.1–3.0 mm for the ICCM and range of 0.1∼4.0 mm for the DCCM; (2) the scale range significantly affects the estimated value of fractal dimension, especially for a Takagi surface; (3) the CCM always provide a higher fractal dimension of the Takagi surface than other three methods do; and the DCCM and RDCCM usually give similar results of the fractal dimensions; (4) the CCM and ICCM may give a result larger than 3.0 for a high fractal dimension surface, while the DCCM and RDCCM can provide a fractal dimension which is located in a range of 2.0–3.0; (5) for a low fractal dimension Takagi surface, the DCCM and RDCCM give accurate results within the ranges determined by small covering scales; for high fractal dimension surfaces, the CCM and ICCM provide very good results at small scales. Acknowledgements The authors would like to express sincere gratitude to all the anonymous reviewers for their comments devoted to improving the quality of our paper. This paper was financially supported by the Major State Basic Research Projects (2011CB201201; 2010CB226802) and National Natural Science Foundation of China (51204113; 51134018). References [1] G.P. Cherepanov, A.S. Balankin, V.S. Ivanova, Fractal fracture mechanics – a review, Eng. Fract. Mech. 51 (1995) 997–1033. [2] O. Mughieda, A.K. Alzo’ubi, Fracture mechanisms of offset rock joints – a laboratory investigation, Geotech. Geol. Eng. 22 (2004) 545–562. [3] T. Belem, F. Homand-Etienne, M. Souley, Fractal analysis of shear joint roughness, Int. J. Rock Mech. Min. Sci. 34 (1997), 130.e131–130.e116. [4] Y. Jiang, B. Li, Y. Tanabashi, Estimating the relation between surface roughness and mechanical properties of rock joints, Int. J. Rock Mech. Min. Sci. 43 (2006) 837–846. [5] W. Kwa´sny, Structure, physical properties and fractal character of surface topography of CVD coatings, J. Mach. Eng. 11 (2011) 134–151. [6] B.B. Mandelbrot, The Fractal Geometry of Nature, Henry Holt and Company, New York, 1983. [7] H. Xie, H. Sun, Y. Ju, Z. Feng, Study on generation of rock fracture surfaces by using fractal interpolation, Int. J. Solids Struct. 38 (2001) 5765–5787. [8] T. Babadagli, K. Develi, Fractal characteristics of rocks fractured under tension, Theor. Appl. Fract. Mech. 39 (2003) 73–88. [9] V. Rasouli, J. Harrison, Assessment of rock fracture surface roughness using Riemannian statistics of linear profiles, Int. J. Rock Mech. Min. Sci. 47 (2010) 940–948. [10] H. Xie, J.-A. Wang, E. Stein, Direct fractal measurement and multifractal properties of fracture surfaces, Phys. Lett. A 242 (1998) 41–50.

T. Ai et al. / Applied Surface Science 314 (2014) 610–621 [11] H.W. Zhou, H. Xie, Direct estimation of the fractal dimensions of a fracture surface of rock, Surf. Rev. Lett. 10 (2003) 751–762. [12] S. Mahovic Poljacek, D. Risovic, K. Furic, M. Gojo, Comparison of fractal and profilometric methods for surface topography characterization, Appl. Surf. Sci. 254 (2008) 3449–3458. ´ S.M. Poljaˇcek, M. Gojo, On correlation between fractal dimension [13] D. Risovic, and profilometric parameters in characterization of surface topographies, Appl. Surf. Sci. 255 (2009) 4283–4288. [14] S. Dhillon, R. Kant, Quantitative roughness characterization and 3D reconstruction of electrode surface using cyclic voltammetry and SEM image, Appl. Surf. Sci. 282 (2013) 105–114. ´ S.M. Poljaˇcek, K. Furic, ´ M. Gojo, Inferring fractal dimension of [15] D. Risovic, rough/porous surfaces – a comparison of SEM image analysis and electrochemical impedance spectroscopy methods, Appl. Surf. Sci. 255 (2008) 3063–3070. [16] C. Ma, W. Wang, S. Li, C. Miao, Q. Zhang, Multifractal, structural, and optical properties of Mn-doped ZnO films, Appl. Surf. Sci. 261 (2012) 231–236. [17] S. Jureˇcka, H. Angermann, H. Kobayashi, M. Takahashi, E. Pinˇcík, Multifractal analysis of textured silicon surfaces, Appl. Surf. Sci. 301 (2014) 46–50. [18] X. Sun, Z. Fu, Z. Wu, Multifractal analysis and scaling range of ZnO AFM images, Phys. A: Stat. Mech. Appl. 311 (2002) 327–338. [19] H.-S. Yu, X. Sun, S.-F. Luo, Y.-R. Wang, Z.-Q. Wu, Multifractal spectra of atomic force microscope images of amorphous electroless Ni–Cu–P alloy, Appl. Surf. Sci. 191 (2002) 123–127.

621

[20] C. Brown, G. Savary, Describing ground surface texture using contact profilometry and fractal analysis, Wear 141 (1991) 211–226. [21] J. Mah, C. Samson, S.D. McKinnon, D. Thibodeau, 3D laser imaging for surface roughness analysis, Int. J. Rock Mech. Min. Sci. 58 (2013) 111–117. [22] H. Xie, J.-A. Wang, M. Kwa´sniewski, Multifractal characterization of rock fracture surfaces, Int. J. Rock Mech. Min. Sci. 36 (1999) 19–27. [23] K.C. Clarke, Computation of the fractal dimension of topographic surfaces using the triangular prism surface area method, Comput. Geosci. 12 (1986) 713–722. [24] W. Kwa´sny, A modification of the method for determination of surface fractal dimension and multifractal analysis, J. Achiev. Mater. Manuf. Eng. 33 (2009) 115–125. [25] Y.H. Zhang, H.W. Zhou, H.P. Xie, Improved cubic covering method for fractal dimensions of a fracture surface of rock, Chin. J. Rock Mech. Eng. 24 (2005) 3192–3196. [26] J. Li, Q. Du, C. Sun, An improved box-counting method for image fractal dimension estimation, Pattern Recognit. 42 (2009) 2460–2469. [27] K. Falconer, Fractal Geometry: Mathematical Foundations and Applications, John Wiley & Sons, Chichester, England, 1990. [28] B. Dubuc, S. Zucker, C. Tricot, J. Quiniou, D. Wehbi, Evaluating the fractal dimension of surfaces, Proc. R. Soc. Lond. A: Math. Phys. Sci. 425 (1989) 113–127. [29] S. Stach, J. Cybo, J. Chmiela, Fracture surface – fractal or multifractal? Mater. Charact. 46 (2001) 163–167.