Computer Estimation of JRC Index Using Function Moments

MENDEL — Soft Computing Journal, Volume 2d, No.gk, .2+2K#2` 2021, Brno, Czech RepublicX ISSN: 1803-3814 (Printed), 2571-3701 (Online) https://doi.org/10.13164/mendel.2021.k.08R Computer Estimation of JRC Index Using Function Moments Dalibor Martiˇsek Institute of Mathematics, Faculty of Mechanical Engineering, Brno University of Technology, Brno, Czech Republic ma[email protected]r.cz Abstract The paper deals with a shape of geological discontinuities. This shape significantly affects the stability of rock massifs. The joint roughness coefficient (JRC) is one of the main shape indicators but the methods of its estimation are based on the empirical analysis of fracture curves in present. We propose a new mathematically correct theory of automatic estimation of the JRC. It shows that determination of the JRC should be based not only on subjective experience, but objective shape characteristics should be used as well. The moment method is one of the possibilities. The principal moments of a fracture surface and the elongation of so called equimomental ellipse can be determined as the possible characteristics of the shape of a fracture surface. The paper introduces a software which is able to reconstruct a 3D profile of a scanned surface and to assign its JRC index automatically. Keywords: Barton Profile, Shape, Fracture Surface, Central Moments, Principal Moments, Equimomental Ellipse, JRC Index. Received: 02 October 2021 Accepted: 07 December 2021 Published: 21 December 2021 1 Introduction The shape of geological discontinuities plays an important role in influencing the stability of rock masses. Many approaches have been used for its determination. The method of Barton and Choubey (see [4]) is well known in geotechnical practice. These authors introduced the method which is able to calculate the shear strength τof rock joints as τ=σnφr+JRC ·log JSC σn where JRC is the joint roughness coefficient, JSC is the joint compressive strength, φris the residual friction angle, and σnis the normal stress. The JRC and JSC indexes are popularly known (see [2,5,6,8,11,12,18,20,22,26,27] for example). However, the JRC index (first proposed in [6]) is the most problematic term in present; the methods of its estimation are based on the empirical analysis of fracture curves – Barton roughness profiles. Barton and Choubey published ten standard joint roughness profiles in a graphical form (see [6]) and suggested a visual comparison of an actual profile with these ten standards. Such assessment, however, can be difficult and very subjective. At present, the profile curves are often described by simple mathematical functions of one variable, for example, a parabola with empirical coefficients (see [23]). However, a rock and its fracture surface is a threedimensional object and should be described in this way. The moment method is very suitable for solving this problem. The so-called complex moments are relatively well known in image processing (see [1], [28], and [25]); the principal geometrical 2D moments are used for an efficient representation of 2D shapes in [24]and[7]. They are suitable for detection of the orientation (see [13], [10], and [29]) and axes symmetry (see [15]) of 2D shapes, and even of potential fields (see [21]). In this paper, we offer a theoretically correct threedimensional description of a fracture surface and its shape. We define the term shape and describe one of several possibilities of its mathematical characterization. A software using this theory for automatic estimation of the JRC index of a scanned rock sample is described as well. 2 Problem Formulation We often discuss the shape of bodies. We say that a body has the shape of a cube, cone, cylinder, or sphere. But what is the shape? It is mostly considered as a qualitative property with a very difficult mathematical characterization, but we can say in general that it is a property invariant to some geometric transformations – translation, rotation, axis symmetry and scaling. In mathematics, there exist a lot of such objects, constructions and expressions with the same properties. One of them are the function moments. We want to describe a three-dimensional property; therefore three-dimensional moments are used. The three-dimensional function moment of the [m;n;p]order is defined as Mf(m;n;p)=Ω xm·yn·zp·f(x;y;z)dxdydz (1) where Ω describes the analyzed object. However, the function moment (1) is not appropriate for a shape characterization because it depends on the location of the coordinate system and its units. But we can define 51 MENDEL — Soft Computing Journal, Volume 2d, No.gk, .2+2K#2` 2021, Brno, Czech RepublicX the center of Ω as the point C=[xC;yC;;zC;]where xC=Mf(1;0;0) Mf(0;0;0) ;C=Mf(0;1;0) Mf(0;0;0) ;zC=Mf(0;0;1) Mf(0;0;0) (2) Now we can define the central moments CMf(m;n;p)=Ωrm(x)·sn(y)·tp(z)·f(r(x); s(y); t(z))dxdydz (3) where r(x)=x−xC s(y)=y−yC t(z)=z−zC (4) With this modification, the central moments are invariant to translation but they are not invariant to other transformations. Therefore, we define the normalized central moments as NCMf(m;n;p)= CMf(m;n;p) [Mf(0; 0; 0)]m+n+p+3 3 (5) The normalized central moments are invariant to translation and scaling, but they are not invariant to rotation. To obtain rotational invariance, we determine the main coordinate system using the eigenvectors of the so-called moment matrix Mf=⎛ ⎝ NCMf(2; 0; 0) NCMf(1; 1; 0) NCMf(1; 0; 1) NCMf(1; 1; 0) NCMf(0; 2; 0) NCMf(0; 1; 1) NCMf(1; 0; 1) NCMf(0; 1; 1) NCMf(0; 0; 2) ⎞ ⎠ i.e., we solve the equation (Mf−λE)h=o(6) We obtain three eigenvectors h1;h2;h3which define the Euler angles α;β;γof the main coordinate axes (the precession, nutation and rotation angle), and three eigenvalues λ1;λ2;λ3asthesizesoftheaxesofthesocalled reference ellipsoid. We rotate the orthonormal base {e1;e2;e3}around the coordinate axes x;y;zby the Euler angles. These rotations are defined as follows: Rα=⎛ ⎝ cos α−sin α0 sin αcos α0 001 ⎞ ⎠ Rβ=⎛ ⎝ cos β0−sin β 01 0 sin β0cosβ ⎞ ⎠ Rγ=⎛ ⎝ cos γ−sin γ0 sin γcos γ0 001 ⎞ ⎠ By transforming the normalized central moments (4) to the main coordinate system we obtain ⎛ ⎝ u(x) v(y) w(z) ⎞ ⎠=Rγ·Rβ·Rα·⎛ ⎝ r(x) s(y) t(z) ⎞ ⎠(7) where r(x);s(y);t(z) are given by (4). Using (5) we can define the principal moments as PMf(m;n;p)= 1 Mf(0;0;0) m+n+p+3 3 · Ω um(x)·vn(y)·wp(z)·f(u(x); v(y); w(z))dxdydz (8) The principal moments are invariant to translation, rotation and scaling, as well as to what we call the shape of a body. Therefore we can say: Two bodies have the same shape if and only if all their principal moments are the same. The principal moments may be therefore used as shape detectors. However, it is necessary to solve two problems: a) A three-dimensional continuous solution outlined above is very complicated. b) Any two fracture surfaces do not have exactly the same shape. In the following text, the problem a) is solved by its discretization and simplification to two dimensions. In the case b), we cannot find exactly the same shape but a “similar” shape only. Therefore, we cannot consider the infinite number of the principal moments but several first only. As is known, the Barton’s JRC standard differentiates ten basic profile curves. For this “shape resolution”, two moments only will be sufficient, represented even with one single number. 3 Materials and Methods All subsequent data used in this paper were acquired by means of a special hardware designed and assembled by prof. Tom´aˇs Ficker from the Faculty of Civil Engineering of our university. All the samples are specimens of limestone (locality Brno-H´ady, Czech Republic). See [16], [17], [19] for more information about its acquisition. Problem discretization: In practice, the data describing analyzed samples are not continuous but discrete sets. Therefore, the integrals in the expressions (1), (3) and (8) can be substituted by sums: Mf(m;n;p)=  [xi;yj;zk]∈Ω xm i·yn j·zp k·f(xi;yj;zk)(9) CMf(m;n;p)=  [xi;yj;zk]∈Ω rm i·sn j·tp k·f(ri;sj;tk) (10) PMf(m;n;p)= 1 [Mf(0; 0; 0)]m+n+p+3 3 ·  [xi;yj;zk]∈Ω um i·vn j·wp k·f(ui;vj;wk) (11) Elimination of the nutation angle: Obviously, geological samples do not have a quite general shape. For example, the part of a sample observed with a microscope or camera can be regarded as a cuboid whose top wall is replaced with the observed profile f(xi;yj). In Fig. 1 we can see an example of an observed block of limestone. Besides the shape that interests us, this block has the “predominant slope”. This slope can be eliminated by a suitable choice of the position of the zaxis, or plane perpendicular to the z-axis – basic plane. 52 MENDEL — Soft Computing Journal, Volume 2d, No.gk, .2+2K#2` 2021, Brno, Czech RepublicX J`iBb2F,g*QKTmi2`g1biBKiBQMgQ7gC_*gAM/2tglbBM;g6mM+iBQMgJQK2Mib Figure 1: Sample of limestone with “predominant left to right slope”. We use the method of least squares for its calculation. We want to interspace a plane z(x;y)=ax +by +c(12) through the measured points [xi;yi;zi], i.e., we have to minimize the function H(a;b;c)= W−1  i=0 H−1  j=0 zij −axi−byj−c2(13) where W×His the resolution of the sample. This problem leads to the system of linear equations a W−1  i=0 x2 i+b W−1  i=0 H−1  j=0 xiyj+c W−1  i=0 xi= W−1  i=0 xi H−1  j=0 zij a W−1  i=0 H−1  j=0 xiyj+b H−1  j=0 y2 j+c H−1  j=0 yj= H−1  j=0 yj W−1  i=0 zij a W−1  i=0 xi+b H−1  j=0 yj+c·W·H= W−1  i=0 H−1  j=0 zij for the unknowns a;b;c. We obtain the plane (12) – for sample from Fig. 1, it is constructed in Fig. 2. We can assume that the z-axis of the profile p(xi;yj)is vertical (see Fig. 3). Our problem is discrete and two dimensional now; expressions (9), (10), (11) are transformed to Mp(m;n)= W−1  i=0 H−1  j=0 xm iyn jp(xi;yj) (14) CMp(m;n)= W−1  i=0 H−1  j=0 rm isn jp(ri;sj) (15) PMp(m;n)= 1 [Mp(0; 0)]m+n+2 2 W−1  i=0 H−1  j=0 um ivn jp(ui;vj) (16) and expressions (2), (4), (6), (7) to xC=Mp(1; 0) Mp(0; 0);yC=Mp(0; 1) Mp(0; 0) (17) ri=xi−xC sj=yj−yC(18)  NCMp(2; 0) −λNCM p(1; 1) NCMp(1; 1) NCMp(0; 2) −λ = 0 (19) ui vj=cos γ−sin γ sin γcos γ·cos α−sin α sin αcos α·ri sj (20) By solving (19), we obtain two eigenvalues λ1;λ2– principal moments PMp(2; 0) ,PMp(0; 2). Principal vectors h1;h2and Euler angles α;γ(β= 0 because of vertical z-axis) will not be needed. 53 MENDEL — Soft Computing Journal, Volume 2d, No.gk, .2+2K#2` 2021, Brno, Czech RepublicX Figure 2: Basic plane detected in sample from Fig. 1. Elimination of the rectangle or square in the top view: 3D data describing a geological sample can be obtained in several ways. Small samples may be scanned with a confocal microscope and transformed to 3D by means of a company software, an accessory of the microscope. In the case of large samples, we can take photographs with a classic or CCD camera and construct the 3D model using methods described in [16], [17], [19], for example. In both cases, the rectangle is in the top view of the reconstruction. It is very unpleasant from the point of view of fracture surface shape identification because the homogenous rectangle or square shape itself has clearly defined principal vectors and moments. The shape of the fracture surface affects these variables very little. To eliminate this unpleasant effect, it is necessary to work with a sample with infinitely many principal axes with the same principal moments in the top view. A circle is such a shape – it is necessary to cut it from the rectangular sample (according to the user selection) – see Fig 4. All the principal functional moments of the profile obtained in this way are given just by the shape of the fracture surface. Moreover, the principal moments of the second order will be sufficient to differentiate the “basic shape” of the fracture surface. Let us sum up the whole algorithm and let us show how to describe a fracture surface by means of just one number. 1. A real sample can be regarded as a cuboid whose top wall is replaced with the profile (xi;yj)obtained by reconstruction of a series of partially focused images. 2. The predominant slope is eliminated using the method of least squares – see (13). 3. The top wall can be considered as smooth now (i.e., horizontal and two-dimensional) and the profile p(xi;yj) can be represented as the density in individual points. 4. A circle is cut from such a surface, according to the user selection. 5. The principal moments of the second order PMp(2; 0) ,PMp(0; 2) are calculated for this circle. The differentness of these two moments can be modeled by means of the so-called equimomental ellipse. It is the ellipse with the same moments PMp(2; 0) ,PMp(0; 2) as our sample. The axes of this ellipse are a=2 PMp(2,0); b=2 PMp(0,2) For our purpose, just one number will be sufficient – the equimomental ellipse elongation, i.e. the number EL =log 2 a b 54 MENDEL — Soft Computing Journal, Volume 2d, No.gk, .2+2K#2` 2021, Brno, Czech RepublicX J`iBb2F,g*QKTmi2`g1biBKiBQMgQ7gC_*gAM/2tglbBM;g6mM+iBQMgJQK2Mib Figure 3: Sample from Fig. 1 with subtracted basic plane from Fig. 2. Figure 4: Circle chosen by user for subsequent calculation of sample shape. Now, for an ideally smooth fracture surface, we obtain a=b, i.e. EL=0. On the other hand, the more the surface of our sample is bumpy (i.e. the more the surface differs from a smooth circle), the grater the elongation is. (Let us note that a=b holds also for an ideally symmetrical surface. However, in practice no real-world sample will be this symmetrical). 4Experiments We acquired 3D reconstruction of two samples of limestone. 3D reconstruction and visualization of these samples can be seen in Fig. 5. We constructed their 3D models using the methods described in [16], [17], [19]. Consequently, the equimomental ellipse elongation of these samples was computed according to the previous section. The random midpoint displacement method represents a de facto standard in natural fractal generation techniques. See [3], [9] for information on this algorithm in 2D, and [14] in 3D. Two of ten of these surfaces can be seen in Fig. 6. We have measured the Equimomental Ellipse Elongation (EEE) of these ten surfaces; results of these measurements are written in Tab. 1. In this way, we obtain a relationship between the interval of the JRC index and the value of the equimomental ellipse elongation in the case of 3D surfaces with Barton profiles. The values of the 3D EEE stated in Tab. 1 correspond to the midpoint of the 2D JRC interval. After that we construct the boundaries of the 3D EEE intervals corresponding to the JRC intervals. 5 Results We have developed a software that works on the principle described above. Let mibethe3DEEEvaluefor the i-th Barton curve stated in Tab. 2. The boundaries of the i-th interval are ai=1 2·(mi−1+mi); i>1 (21) (lower boundary) and bi=1 2·(mi+mi+1); i<10 (22) (upper boundary). We put a0= 0 and b10 →∞. The program provides the resulting value of the JRC of the chosen area (for 55 MENDEL — Soft Computing Journal, Volume 2d, No.gk, .2+2K#2` 2021, Brno, Czech RepublicX Figure 5: Samples A, B in which elongations were measured. Figure 6: Computer modeling of surface with JRC = 0-2 and 18-20 (random midpoint displacement method). Table 1: Intervals of JRC indexes and corresponding intervals of 3D EEE. IJRC 3D EEE interval – see (21), (22) aimibi 10-2 0.000 0.072 0.080 22-4 0.080 0.088 0.098 34-6 0.098 0.107 0.118 46-8 0.118 0.128 0.142 58-10 0.142 0.155 0.168 610-12 0.168 0.181 0.195 712-14 0.195 0.208 0.220 814-16 0.220 0.232 0.245 916-18 0.245 0.258 0.266 10 18-20 0.266 0.273 →∞ the samples in Fig. 1 we obtain the elongation 0.109, i.e. the JRC index is 6. 6 Conclusion Visual comparison of an observed profile with the Barton standards can be difficult, very subjective and always uncertain. The method presented above is objective and reliable in the case of quality input data. The elongations for different samples are different and the difference between the values for the individual samples can be taken as a measure of the difference between jaggednesses of their surfaces (a small difference between samples A and B indicates “almost the same shape”) and between shear strengths too. There exists a relationship between the equimomental ellipse elongation and the joint roughness coefficient. This relationship may be used for the automatic software estimation of the JRC index of various rock surfaces. This fact was demonstrated by means of a functional software based on this principle. Our software is able to process speciments with different sizes. It means, it is able to provide 3D reconstruction and principialy also JRC estimation of the speciment which size is several tens of micrometers only. However, JRC estimation is reliable in the case of the speciments which size is at least several centimeters. 3D reconstruction of smaller samples may be used for other purposes (morphological analysis of fracture surfaces of steel or building materials for example). The interval of the JRC is assigned to the interval ofthe3DEEEvalueinthisway. Thisassignment can be seen in Tab. 2. This table enables to assign the JRC index to each calculated equimomental ellipse elongation value. The reliability of our method can be decreased by additional noise in the input data; a large share of additive noise (low signal-to-noise ratio) can increase the 56 MENDEL — Soft Computing Journal, Volume 2d, No.gk, .2+2K#2` 2021, Brno, Czech RepublicX J`iBb2F,g*QKTmi2`g1biBKiBQMgQ7gC_*gAM/2tglbBM;g6mM+iBQMgJQK2Mib Table 2: Barton roughness 2D profiles (original on the left – see [2], processed by Image processing method on the right) and equimomental ellipse elongation of corresponding 3D profiles generated by means of software random midpoint displacement method. Curve Typical 2D Roughness Profile 2D JRC 3D EEE 10-2 0.072 2 2-4 0.088 3 4-6 0.107 4 6-8 0.128 5 8-10 0.155 6 10-12 0.181 7 12-14 0.208 8 14-16 0.232 9 16-18 0.258 10 18-20 0.273 estimation of the JRC and reduce the reliability of this estimation. Acknowledgement: The author acknowledges support from Private Institute of Applied Mathematics, Slapanice, Czech Republic. The author would like to thank prof. Tom´aˇs Ficker from the Faculty of Civil Engineering, and ass. prof. Pavel ˇ Starha from the Faculty of Mechanical Engineering (both from Brno University of Technology) for the provided data. References [1] Abu-Mostafa, Y. S., and Psaltis, D. Image normalization by complex moments. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1 (1985), 46–55. [2] Amanloo, F., and Hosseinitoudeshki, V. The effect of joint roughness coefficient (jrc) and joint compressive strength (jcs) on the displacement of tunnel. International Research Journal of Applied and Basic Sciences, Science Explorer Publications 4, 8 (2013), 2216–2224. [3] Barnsley, M. F., Devaney, R. L., Mandelbrot, B. B., Peitgen, H.-O., Saupe, D., Voss, R. F., Fisher, Y., and McGuire, M. The science of fractal images, vol. 1. Springer, 1988. [4] Barton, N. Review of a new shear-strength criterion for rock joints. Engineering Geology 7,4 (1973), 287–332. [5] Barton, N. Shear strength criteria for rock, rock joints, rockfill and rock masses: Problems and some solutions. Journal of Rock Mechanics and Geotechnical Engineering 5, 4 (2013), 249–261. [6] Barton, N., and Choubey, V. D. The shear strength of rock joints in theory and practice. Rock mechanics 10 (1977), 1–54. [7] Crespo, J. F., Crespo, J. F., Lopes, G. A., and Aguiar, P. M. Principal moments for efficient representation of 2d shape. In 2009 16th IEEE International Conference on Image Processing (ICIP) (2009), IEEE, pp. 1085–1088. [8] Du, S., Hu, Y., Hu, X., and Guo, X. Comparison between empirical estimation by jrc-jcs model and direct shear test for joint shear strength. Journal of Earth Science 22, 3 (2011), 411–420. [9] Fournier, A., Fussell, D., and Carpenter, L. Computer rendering of stochastic models. Communications of the ACM 25, 6 (1982), 371– 384. [10] Ha, V. H., and Moura, J. M. Efficient 2d shape orientation. In Proceedings 2003 International Conference on Image Processing (Cat. No. 03CH37429) (2003), vol. 1, IEEE, pp. I–225. [11] Han, F.-s., and Tang, C.-a. Numerical investigation for anisotropy of compressive strength of rock mass with multiple natural joints. Journal of Coal Science and Engineering (China) 16,3 (2010), 246–248. [12] Li, Y., Wang, J., Jung, W., and Ghassemi, A. Mechanical properties of intact rock and fractures in welded tuff from newberry volcano. In Proceedings of 37th Workshop on Geothermal Reservoir Engineering, Stanford, CA (2012), vol. 30. [13] Lin, J.-C. Universal principal axes: an easy-toconstruct tool useful in defining shape orientations for almost every kind of shape. Pattern Recognition 26, 4 (1993), 485–493. [14] Mandelbrot, B. B., and Mandelbrot, B. B. The fractal geometry of nature, vol. 1. WH freeman New York, 1982. [15] Marola, G. On the detection of the axes of symmetry of symmetric and almost symmetric planar images. IEEE Transactions on Pattern Analysis and Machine Intelligence 11, 1 (1989), 104–108. [16] Martiˇ sek, D. The two-dimensional and threedimensional processing of images provided by conventional microscopes. Scanning: The Journal of Scanning Microscopies 24, 6 (2002), 284–296. [17] Martiˇ sek, D., and Druckm¨ ullerov´ a, H. Multifocal image processing. Mathematics for Applications 3, 1 (2014), 77–90. [18] Martisek, D., and Prochazkova, J. The analysis of rock surface asperities. MENDEL Journal 24, 1 (2018), 135–142. 57 MENDEL — Soft Computing Journal, Volume 2d, No.gk, .2+2K#2` 2021, Brno, Czech RepublicX [19] Martiˇ sek, D., Proch´ azkov´ a, J., and Ficker, T. High-quality three-dimensional reconstruction and noise reduction of multifocal images from oversized samples. Journal of Electronic Imaging 24, 5 (2015), 053029. [20] Nakagawa, M., Jiang, Y., Kawakita, M., Yamada, Y., and Akiyama, Y. Evaluation of mechanical properties of natural rock joints for discontinuous numerical analysis. In Proc. ISRM Int. Symp. 3rd ARMS (2004), Millpress. [21] Prasad, V. S. N., and Yegnanarayana, B. Finding axes of symmetry from potential fields. IEEE Transactions on Image Processing 13,12 (2004), 1559–1566. [22] Prudencio, M., and Van Sint Jan, M. Strength and failure modes of rock mass models with non-persistent joints. International Journal of Rock Mechanics and Mining Sciences 44,6 (2007), 890–902. [23] Rafek, A. G., and Goh, T. L. Correlation of joint roughness coefficient (jrc) and peak friction angles of discontinuities of malaysian schists. Earth Science Research 1, 1 (2012), 57. [24] Rodrigues, J. J., Aguiar, P. M., and Xavier, J. M. Ansig—an analytic signature for permutation-invariant two-dimensional shape representation. In 2008 IEEE Conference on Computer Vision and Pattern Recognition (2008), IEEE, pp. 1–8. [25] Shen,D.,Ip,H.,Cheung,K.,andTeoh, E. K. Symmetry detection by generalized complex (gc) moments: a close-form solution. IEEE Transactions on Pattern Analysis and Machine Intelligence 21, 5 (1999), 466–476. [26] Singh, M., and Singh, B. High lateral strain ratio in jointed rock masses. Engineering Geology 98, 3-4 (2008), 75–85. [27] Starha, P., Prochazkova, J., and Martisek, D. The reconstruction of the object surface using confocal microscope with hyperchromatic lens. MENDEL Journal 24, 1 (2018), 129– 134. [28] Teh, C.-H., and Chin, R. T. On image analysis by the methods of moments. IEEE Transactions on pattern analysis and machine intelligence 10,4 (1988), 496–513. [29] ˇ Zuni´ c, J., Kopanja, L., and Fieldsend, J. E. Notes on shape orientation where the standard method does not work. Pattern Recognition 39, 5 (2006), 856–865. 58