Arch. Min. Sci., Vol. 57 (2012), No 4, p. 921932 Electronic version (in color) of this paper is available: http://mining.archives.pl DOI 10.2478/v10267-012-0061-y MASOUD SOLEYMANI SHISHVAN*, JAVAD SATTARVAND** MODELING OF ACCURATE VARIABLE SLOPE ANGLES IN OPEN-PIT MINE DESIGN USING SPLINE INTERPOLATION MODELOWANIE ZMIENNEGO KTA NACHYLENIA STOKU W PROJEKTOWANIU KOPALNI ODKRYWKOWYCH ZA POMOC INTERPOLACJI FUNKCJAMI SKLEJAJCYMI (METOD SPLINE'ÓW) In this paper a new method of modeling variable slope angles has been presented based on the spline interpolation method. Slope angle modeling and defining precedency of the blocks are the vital parts of almost any open pit optimization algorithm. Traditionally heuristic patterns such as 1:5 or 1:9 have been used to generate slope angles. Cone template based models were later employed in developing variable slope angles. They normally use a linear interpolation process for determination of slope angles between the given directions which leads to sharp and non-realistic pits. The other elliptical alternatives suffer from having limitations in defining slope angles in non-geographical directions. The method is capable to consider any number of slope angles in any desired direction as well as creating quite accurate and realistic pit shapes. Three major types of the spline interpolation including cubic, quadratic and cardinal are tested, however, the cubic form is preferred due to more realistic outcomes. Main steps of the method are described through a numerical case study. Keywords: open pit optimization, variable slope angles, spline interpolation, production planning W pracy zaprezentowano now metod modelowania zmiennego kta nachylenia gruntu w oparciu o metod interpolacji funkcjami sklejajcymi (metoda spline'ów). Modelowanie kta nachylenia stoku i prognozowanie kolejnoci wybierania to kluczowe elementy algorytmu optymalizacyjnego. Tradycyjne modele heurystyczne oparte o wzorce 1:5 lub 1:9 wykorzystane zostaly do wygenerowania któw nachylenia stoku. Do wygenerowania zmiennych któw nachylenia wykorzystano modele stokowe. Procedura taka zasadniczo zaklada wykorzystanie interpolacji liniowej dla okrelenia kta nachylenia pomidzy dwoma kierunkami, co prowadzi moe do zaprojektowania bardzo stromych i nierealistycznych ksztaltów odkrywek. Alternatywne rozwizania, wykorzystujce modele eliptyczne, maj inne ograniczenia mianowicie okrelaj one kty nachylenia w kierunkach innych ni geograficzne. Za * ** MINING ENGINEERING, SAHAND UNIVERSITY OF TECHNOLOGY, TABRIZ, IRAN, E-mail: M_soleymani@sut.ac.ir, TEL: 0098 412 3444312 SAHAND UNIVERSITY OF TECHNOLOGY, TABRIZ, IRAN, E-mail: email@example.com, TEL: 0098 412 3459267 ACADEMIC MEMBER, INSTITUTE OF SURFACE MINING AND DRILLING, RWTH AACHEN UNIVERSITY, GERMANY, EMAIL: SATTARVAND@BBK3.RWTH-AACHEN.DE, TEL: 0049 241 8095683 pomoc tej metody uwzgldni mona dowoln liczb któw nachylenia w dowolnym kierunku a take wygenerowa dokladne i realistyczne ksztalty odkrywek. Przetestowano trzy procedury interpolacyjne: z zastosowaniem funkcji szeciennych, kwadratowych i kardynalnych. Zdecydowanie najkorzystniejsze i najbardziej realistyczne wyniki uzyskuje si przy zastosowaniu funkcji szeciennych. Glówne etapy stosowanej metody wyjanione zostaly przy pomocy przykladu numerycznego. Slowa kluczowe: optymalizacja odkrywki, zmienne kty nachylenia stoku, interpolacja funkcjami sklejajcymi (metod spline'ów), planowanie wydobycia 1. Introduction The extensive application of computers in the field of production planning of open-pit mines has engaged researchers to develop superior and comprehensible algorithms for answering the problem as fast as possible and requiring less computer and human resources. However even after about 40 years, the field still needs to develop more powerful methods to define Ultimate Pit Limit (UPL), optimum cut-off grade, optimum mine schedule, mine life and capacities. They also need to undertake all risks and uncertainties involved in estimation of geological grades and market condition. In early studies (Lemieux, 1979; Pana, 1965; Williams, 1974) developed the moving cone algorithm for designing the outline of the final pit shape. Several studies were carried out later to combine 2D sections into 3D pit shapes (Johnson & Sharp, 1971; Wright, 1987). (Koenigsberg, 1982) and (Wilke & Wright, 1984) succeeded in directly applying dynamic programming to solve a 3D pit design problem. Lerchs and Grossmann's graph theory based algorithm (Lerchs & Grossmann, 1965) might be one of the most well-known and utilized algorithms in the field of open-pit optimization in order to formulate UPL problem. Attempts continued afterwards to develop more efficient algorithms for this problem (Hochbaum, 2001; Huttagosol & Cameron, 1992; Yegulalp & Arias, 1992; Zhao & Kim, 1992). Subsequent studies focused on more general problem rather than the UPL which was the production planning problem. This challenging problem tries to determine the excavation time and the destination of each exploited block. Early researches attempted to solve the mathematical model of the problem (Caccetta, Kelsey, & Giannini, 1998; Dagdelen & Johnson, 1986; Ramazan, Dagdelen, & Johnson, 2005). Recently some studies have been reported on using metaheuristic algorithms too. (Denby & Schofield, 1996) tried to use genetic algorithm, (Kumral & Dowd, 2005) proposed another metaheuristic algorithm based on simulated annealing and newly (Sattarvand, 2009) presented a new algorithm based on ant colony optimization. Lately (Sayadi, Fathianpour, & Mousavi, 2011) used a new artificial neural network method for optimizing open pit production planning. The slope angle of open pit mine is one of the most essential features during whole mine planning and design process. It strongly affects the final shape of the mine and the best layout of the access roads, waste and low grade ore dumps, stockpiles, processing plant, and other surface facilities. Consequently the minable reserves and the amount of ore and waste to be removed during mine life are highly dependent to the slope angles; and almost all of the open pit production planning algorithms are somehow connected to the definition of the slope angles. For example the well-known Lerchs and Grossmann's graph theory based algorithm uses a list of blocks which have to be removed prior to mining of a certain block. This list is normally determined according to the proposed slope angles. Similarly, most of the mathematical scheduling algorithms have a sequencing constraint in their optimization modeling structure that is also defined based on the slope angles. One of the characteristic features of the open pit mining is that the material, which is billions of years old, has been affected by different pressures, displacements, temperatures, chemical processes and tectonic forces cannot be expected to show homogenous stability behavior in all directions. In other words, every small part of the deposit and surrounding rocks would behave quite differently. Authors distinguish two types of the variation in slope angles' pattern of a deposit. The first form, called variable slope angles, considers a set of different slope angles in several horizontal directions. This pattern supposed to be similar for all blocks of the model. Whereas in second type, called multiple variable slope angles, the slope pattern is considered to be individually different for each block of the model based on the rock type. In this paper a new interpolation method has been proposed for accurately modeling of the variable slope angles in open pit planning and design. 2. Review of slope angle modeling approaches Early modeling routines were based on the special block configurations such as 1:5 or 1:9 patterns, i.e., any given block is considered accessible if 5 or 9 blocks on top have been removed before. These approaches were suffering from the creation of higher and lower slope angles than desired. Later (Lipkewich, 1969) proposed a knight move pattern to estimate the conical extensions on the surface by proposing to use a 1:5:9 pattern which was able to create enhanced slopes. However the major disadvantage of all routines was that the created slope angle was dependent to the block dimensions. The problem has been partly overcome by introduction of the idea of cone template. Nowadays cone template is the most frequently used routine in modeling of the variable slope angles. It could be defined as a cone that is constructed above a certain mining block, containing all blocks which have to be removed before extraction of that block. (Chen, 1976) attempted to insert variable slope angles concept within Lerchs and Grossmann algorithm by using this concept. (Dowd & Onur, 1993) utilized the cone template concept in their proposed algorithm. Whittle has repeatedly reported the incorporation of the variable pit slopes by linear interpolation, (Alford, Whittle, & 1986). The basic concept of incorporation of the variable slopes into Lerchs and Grossmann algorithm has been addressed plainly by (Khalokakaie, Dowd, & Fowell, 2000). They assumed that the slope angles are defined in four principal directions (north, south, east and west) and then using an elliptic equation for connecting cone corners on each level a smooth cone template is generated. However, in real cases it is almost unpractical to use only four slope angles in four principal directions. It gets even more unreliable when the orientation of the block model does not coincide the geographical principal directions. Current paper describes a new variable slope angles modeling method using the spline interpolation theory. This method is based on the conventional cone template technique and has the capability to consider unlimited number of variable slope angles in different directions. It defines the cone template by outlining its horizontal sections in each level of block model as a set of closed spline curves. Each spline curve attempts to connect the corners of the slope lines on a level with the best fitted spline function, Figure 1b. Following, after description of the basic theory of polynomial interpolation and spline function, the results of the application of the proposed method on a case study has been presented. a) Fig. 1. Cone template constructed by (a) linear and (b) spline interpolation b) 3. Basics As described, the modeling of variable slope angles of open pits could be expressed as defining a smooth 2D curve over intersection points of the slope lines in different azimuths. These points are employed as data points to construct an interpolation curve by passing the curve through them, (Salomon, 2006). 3.1. Parametric Curves Ordinarily curves are represented implicitly in form of F(x, y) = 0. In practical applications where the complex curves are needed and the function is unknown, a curve would be expressed parametrically as P(t) = (f(t), g(t)). In this representation, the functions f and g are the x and y coordinates of any point on the curve depending on the variable t that varies over a certain interval, normally [0, 1]. The slope of a two-dimensional parametric curve could be expressed as Pyt (t ) dy dy dx ) = t , (Salomon, 2006). =( dx dt dt Px (t ) 3.2. Polynomial interpolation A polynomial of degree n in x is the function: Pn (x) = i= 0 å ai x i = a0 + a1 x + a2 x 2 + ... + an x n (1) where ai are the coefficients of the polynomial. The function has n + 1 coefficients. The cubic polynomial is the simplest curve that would fit complex shapes and has the form of P3 (t) = At 3 + Bt 2 + Ct + D. The function would be determined by finding its four unknown coefficients using four equations based on the four known quantities denoted here by G1 through G4. Therefore a cubic polynomial segment could be expressed as the product: æ m11 ç m21 P (t ) = (t 3, t 2, t ,1) ç ç m31 ç è m41 m12 m22 m32 m42 m13 m 23 m 33 m43 m14 æ æ G1 æ çç ç m 24 G2 ç ç ç = T (t ) × M × G m34 ç ç G3 ç çç ç m44 è è G4 è (2) where M is the basis matrix that depends on the employed method and G is the geometry vector, consisting of the four given quantities (Salomon, 2006). 3.3. Hermite curve Hermite interpolation is a polynomial interpolation that is based on two points and two tangent vectors. In other words, it is a parametric curve that computes the curve segment which starts at P1, going in direction P1t and ends at P2 moving in direction P2t. For a cubic polynomial the function and tangent vectors could be algebraically represented as: T P( t) = at 3 + bt 2 + ct + d = ( t 3, t 2, t,1)(a , b , c, d ) = T ( t) A (3) (4) Pt (t) = 3at 2 + 2bt + c Now the known conditions could be used to construct equations for defining the four unknown coefficients a, b, c and d as following: P(0) = P1 P(1) = P2 Pt (0) = P1t Pt (1) = P2t The solutions are easily obtained as: (5) (6) (7) (8) a = 2 P1 - 2 P2 + P1t + P2t , b = - 3P + 3P2 - 2P1t - 2P2t , 1 c = P1t , d = P1 Substituting these values into Equation (3) gives: t t P (t ) = ( 2P1 - 2P2 + P1t + P2 ) t 3 + ( -3P1 + 3P2 - 2P1t - 2 P2 ) t 2 + P1t t + P1 (9) (10) Which could be rearranged in matrix notation as: 1ææ P æ 1 æ 2 -2 1 ç çç ç - 3 3 - 2 -1 ç ç P2 ç P ( t ) = T .H .B = (t 3, t 2, t ,1) ç ç0 0 1 0 ç ç P1t ç ç çç ç 0 0 0 è ç P2t ç è1 è è (11) where the matrix H is called the Hermite basis matrix, (Salomon, 2006). Hermite basis matrices could be found similarly for higher degree polynomials. 4. Spline Interpolation A spline is a set of polynomials of degree k (Hermite curves) that are smoothly connected at certain data points. Smoothly connection of the curve segments at each data point implies that their tangent vectors (first derivatives) to be equal on these points, (Salomon, 2006). Having n data points, numbered from P1 to Pn, there will be n 1 segments in the model; however, a closed spline has an extra curve segment from Pn to P1 that closes the curve. It is usually convenient to define two additional points Pn + 1 = P1 and Pn + 2 = P2 to define a closed spline. The degree of polynomial (k) is normally considered as 2, 3 or 5 corresponding to quadratic, cubic and quantic splines respectively. Other degrees are not recommended due to adding unnecessary complexity to the calculations. The cubic spline k = 3 that was originally introduced by (Ferguson, 1964) ends to the equations system as: é1 ê0 ê0 ê1 ê4 ë 0 ù 0 ú ú 1 ú 4 ú 1 û (n ´ n) é P1t ù ê ú é 3( P - P ) ù 3 1 ê P2t ú ê ú ê t ú ê 3( P4 - P2 ) ú P3 ú ê ê ú = ê ú ê ú ê ú ê3(Pn +1 - Pn -1 )ú ê Pnt -1 ú ê ú ê ú ë 3(Pn + 2 - Pn ) û t ê Pn ú ë û (12) By separately substitution of x and y coordinates in P1 to Pn variables and solution of the system, the values of P1t to Pnt could be found and subsequently inserted in formula (10) to reach following equation, (Salomon, 2006): P ( t ) = Pk + Pkt × t + (3 (Pk +1 - Pk ) - 2Pk - Pk +1) × t + (2 (Pk - Pk +1) + Pk + Pk +1) × t k t t 2 t t 3 (13) x and y coordinates of the spline curve points could be determined by giving different values to the t from distance of [0, 1]. For a quadratic spline curve, the geometrically representation of segments are as (Salomon, 2006): Pk ( t) = (Pk +1 - Pk - Pkt ) × t 2 + Pkt × t + Pk (14) Whose required parameters could be found by solution of following matrix similarly. é1 ê0 ê0 ê0 ê1 ë ù 0 ú ú 0 ú 1 ú 1 û (n ´ n) é P1t ù ê t ú ê P2 ú ê t ú ê P3 ú = ê ú ê ú ê Pnt-1 ú ê t ú ê Pn ú ë û é ê ê ê ê ê ê ë ù ú 2 ( P3 - P2 ) ú ú ú 2 (Pn - P -1 ) ú n 2 (Pn +1 - Pn ) ú û 2 ( P2 - P1) (15) The cardinal spline is another popular spline interpolation which normally supplies a series of local regulators to control the tension of the curve by modifying the magnitudes of the tangent vectors. The cardinal spline of n given points could be considered as a series of segments, each is defined based on four points. Therefore to construct a cardinal spline, the set of points should be organized into n - 3 overlapping groups of four consecutive points as following: [P, P2 , P3, P ] , [P2 , P , P4, P ] , [P , P , P , P6 ] ,...,[P - 3, P -2 , P -1, P ] ,...,[Pn , P +1, P +2 , P +3] 1 4 3 5 3 4 5 n n n n n n n (16) Hermite interpolation is then applied to construct a curve segment P(t ) for each group as: æ 2 -2 ç -3 3 3 2 P ( t ) = ( t , t , t ,1) ç ç0 0 ç 0 è1 æ -s ç 2s = ( t 3, t 2, t,1) ç ç -s ç è 0 2-s s -3 0 1 1 -2 1 0 1æ ç -1 ç 0ç ç 0è s-2 æ P2 æ ç ç P3 ç ç ç s(P3 - P ) ç 1 ç ç è s(P4 - P2) è -s s æ ç ç 0ç ç 0è æPæ 1 ç ç P2 ç ç ç P3 ç ç ç è P4 è (17) 3 - 2s s 0 Where s is a parameter that controls the magnitude of the spline tension. Spline curve corresponding to s = is called the Catmull Rom or zero tension spline (Catmull & Rom, 1974). Decreasing s from 1/2 to 0 leads to reduction in the magnitude of the tangent vectors and resulting straight segments. In contrast increasing s from 1/2 to 1 results in a curve with more slack at the data points (Salomon, 1999). 5. Implementation of the spline function Major steps of an accurate cone template construction procedure for a given set of azimuths and slope angles are demonstrated in Figure 2. In this process, the shape of a cone template is defined by determining its horizontal intersections (Figure 1). Figure 2 shows that each intersection of the cone template on upper levels of block model is defined by a splines curve. Construction of each spline is accomplished through the following stages: · For each slope direction an imaginary slope line is considered and the coordinates of its intersection points on each level of the block model are determined by simple trigonometric relations. These points are thus called corner points. · Spline interpolation matrices are generated and solved for the set of corner points on each level. · Giving values from zero to one to the variable t in the obtaining geometrical functions, coordinates of a series of points on each spline segment could be calculated. They are called border points in this paper. This is done separately for x and y coordinates. · Border points are sorted according to their azimuth from 0° to 360°. Then all blocks of the desired level are checked to be inside or outside of the cone template by comparing their distance from the cone apex to that of the border point in their direction. procedure Determining of pit shape input azimuths and slope angles for i = 1 to level number do Find corner points Create spline matrix and solve Calculate the border points Find in-cone blocks End-for end-procedure Fig. 2. Process of creating splines A computer program is developed in visual studio (C#) 2010 express edition programming environment for implementation of the calculations. Results of the application of the algorithm with cubic, quadratic and cardinal spline (with s = 1) on six given corner points are shown in Figure 3. Examinations show that the results of quadratic and cardinal splines do not reveal realistic in most of the cases, see Figure 3. Nevertheless, cubic spline is, indeed, able to produce perfectly smooth and representative curves between points. a) b) c) Fig. 3. Example results of a) quadratic spline b) cardinal spline c) cubic spline 6. Case study In order to test the process, a hypothetical block model containing 31×31×9 blocks with block dimension of 10×10×10 m parallel to the principle geographical directions is considered. Slope configuration is also considered to be as Table 1. A sample block is considered in the middle of the lowest level of the block model and accurate cone template is generated using cubic spline interpolation. Initially the set of corner points are defined on each level and spline matrix is constructed and solved for each set. For example the coordinates of the corner points on first level could be expressed as Table 2. The matrices of spline interpolation for these corner points are as following. TABLE 1 Input data set Point 1 Azimuth Slope angle Point 2 Point 3 Point 4 Point 5 TABLE 2 Corner points on the first level of the case study Point_Number P1 P2 P3 P4 P5 X value (m) Y value (m) 10.7090 0.5612 8.1601 6.3754 6.5982 9.4233 11.7365 2.0695 For x coordinates, equation (12) becomes: é1 ê0 ê ê0 ê ê1 ê4 ë é t ù 4 1 0 0 ù ê P1 (x) ú é 3(P3 (x) - P1 (x))ù 1 4 1 0 ú ê P2t (x) ú ê 3(P4 (x) - P2 (x))ú úê ú ú ê 0 1 4 1 ú ê P t (x) ú = ê 3(P5 (x) - P3 (x))ú 3 ú ê ú 0 0 1 4 ú ê P t (x) ú ê 3(P1 (x) - P4 (x))ú 4 ê ú ú ê ú 1 0 0 1 û ê t ú ë 3(P2 (x) - P5 (x))û P5 (x) û ë (18) this equation is solved with the values of Table 2, thus we have: é P1t ( x) ù 17.9385 ù ê t ú é ê ú ê P2 ( x) ú -0.0855 ú ê t ú ê P3 ( x) ú = ê 0.4247 ú ê ê t ú ê -13.9456ú ú ê P4 ( x) ú ê ê P t ( x) ú ê -4.3322 ú ë û ë 5 û Similar equation is written for y coordinate and its solution becomes: (19) é P1t ( y) ù -2.9836 ù ê t ú é ê ú ê P2 ( y) ú -10.3692 ú ê t ú ê P3 ( y) ú = ê -5.0527 ú ê ê t ú ê 3.9939 ú ú ê P4 ( y) ú ê ê P t ( y) ú ê 14.4117 ú ë û ë 5 û (20) Substitution of the solution in the geometrical formula (i.e. equation (13)) of the spline results in: ì P ( x ) ( t) = P ( x) + P t ( x ) × t + (3 (P ( x ) - P ( x)) - 2 P t (x) - P t ( x )) × t 2 k k k k +1 k k k +1 ï ï í î + (2 (Pk ( x) - Pk +1 ( x )) + Pkt (x) + Pkt+1 (x)) × t 3 Pk ( y) (t ) = Pk ( y) + Pkt ( y) × t + (3 ( P +1 ( y) - Pk ( y)) - 2Pkt ( y) - P t+1 ( y)) × t 2 k k + ( 2(P ( y) - Pk +1 ( y)) + Pkt ( y) + Pkt+1( y)) × t 3 k (21) Giving different values to variable t from the distance of [0,1], border points of the spline are determined and converted to the polar coordinate system to produce the distance of border points from center block. Selected numbers of the results are indicated in Table 3 for azimuth distances of 10 degrees. Final result is shown in Figure 4. TABLE 3 Distance of the border points from center block in different azimuths Azimuths Radius Azimuths Radius Azimuths Radius Azimuths Radius Fig. 4. Final results of implementation of spline interpolation on case study 7. Conclusion Presented algorithm is potentially able to enhance the quality of generated pits in almost the entire of open pit production planning procedures. Customarily linear regression method is used in current commercial software for finding of the unknown slope angles between the given directions that leads to non-realistic polygonal shapes of the pit. In fact, the proposed procedure substitutes the sharp polygonal intersections of the cone templates by a set of smooth spline curves on the block model levels. It can consider any number of slope angles in any direction. Considering the fact that the cone template construction is only done one time for each slope region of the model during whole open pit optimization process, increasing the calculation time due to the application of the presented algorithm will be negligible (only some milliseconds). The process is also capable to be utilized in all open pit optimization algorithms for designing ultimate pit limits or production scheduling.
Archives of Mining Sciences – de Gruyter
Published: Dec 1, 2012