Arch. Min. Sci., Vol. 61 (2016), No 3, p. 635649 Electronic version (in color) of this paper is available: http://mining.archives.pl DOI 10.1515/amsc-2016-0045 SAEED SOLTANI-MOHAMMADI*, MOHAMMAD SAFA* SYMULACJA ALGORYTMU OPTYMALIZACYJNEGO PROCESU ODPRANIA DLA AUTOMATYCZNEGO DOPASOWANIA MODELU WARIOGRAMU Fitting a theoretical model to an experimental variogram is an important issue in geostatistical studies because if the variogram model parameters are tainted with uncertainty, the latter will spread in the results of estimations and simulations. Although the most popular fitting method is fitting by eye, in some cases use is made of the automatic fitting method on the basis of putting together the geostatistical principles and optimization techniques to: 1) provide a basic model to improve fitting by eye, 2) fit a model to a large number of experimental variograms in a short time, and 3) incorporate the variogram related uncertainty in the model fitting. Effort has been made in this paper to improve the quality of the fitted model by improving the popular objective function (weighted least squares) in the automatic fitting. Also, since the variogram model function (£) and number of structures (m) too affect the model quality, a program has been provided in the MATLAB software that can present optimum nested variogram models using the simulated annealing method. Finally, to select the most desirable model from among the single/multi-structured fitted models, use has been made of the cross-validation method, and the best model has been introduced to the user as the output. In order to check the capability of the proposed objective function and the procedure, 3 case studies have been presented. Keywords: Automatic variogram fitting, Geostatistics, Optimization, simulated annealing Dopasowanie modelu teoretycznego do eksperymentalnego wariogramu jest kluczowym zagadnieniem w badaniach geostatystycznych poniewa jeli parametry modelu wariogramu obarczone s niepewnoci, to otrzymamy znaczny rozrzut wyników oblicze i symulacji. Pomimo, e najpopularniejsz metoda dopasowania jest dopasowanie `na oko', w niektórych przypadkach wykorzystuje si automatyczne metody dopasowania modelu oparte na zasadach geostatystyki i optymalizacji w celu: 1) dostarczenia podstawowego modelu do dopasowania `na oko'; 2) dopasowania modelu do wikszej iloci eksperymentalnych wariogramów w krótkim okresie czasu; 3) uwzgldnienia niepewnoci zwizanej z wariogramem w dopasowaniu modelu. W pracy podjto prób poprawy jakoci dopasowania modelu poprzez wprowadzenie zmodyfikowanej popularnej funkcji celu (waone najmniejsze kwadraty) do au- UNIVERSITY OF KASHAN, DEPARTMENT OF MINING ENGINEERING, KASHAN, IRAN tomatycznego dopasowania. Ponadto, poniewa funkcja modelu wariogramu (L) i ilo struktur (m) ma take wplyw na jako modelu, opracowano program w rodowisku MATLAB który podaje optymalne modele wariogramu w oparciu o metod symulacji odprania. W czci kocowej wybrano najkorzystniejszy model sporód modeli dopasowania z wykorzystaniem metody walidacji krzyowej i najlepszy model przedstawiany jest uytkownikowi. W celu zbadania moliwoci stosowania proponowanej funkcji celu i przedstawionej procedury, zaprezentowano trzy studia przypadku. Slowa kluczowe: automatyczne dopasowanie wariogramu, geostatystyka, optymalizacja, symulacje procesu odprania 1. Introduction A key stage in geostatistical studies is the structural analysis of the regionalized variable to describe the spatial variations of the desired variable (grade, thickness, accumulation, etc.) (Chiles & Delfiner, 2012). It is the variogram that characterizes the spatial variability of a regionalized variable (Pardo-Igúzquiza, 1999). Generally, the variogram of a regionalized variable is unknown and the experimental variogram should be calculated using the data gathered from the regionalized variable. Using the moments method, the experimental variogram is calculated as follows(Matheron, 1965): ^ h 1 2N h N h z xi k 1 z xj (1) ^ where (h) is the value of the experimental variogram for the vector h= |xi xj|, N(h) is the number of data pairs for h, and z(xi) and z(xj). are the observed or experimental value at points xi and xj. After calculating the experimental variogram, it is necessary that an appropriate model be fitted to it; this is a basic step in geostatistics (Journel & Huijbregts, 1978, Christakos 1984) which largely affects the computation of the krigging weights (Gorsich & Genton, 2000) and, hence, it affects the efficacy of the maps obtained from the krigging estimator and the geostatistical simulation (Genton, 1998a). Although fitting by eye (a procedure in the graphical geostatistical software based on trial and error) is still the most common method of fitting a model to an experimental variogram, with the increasing use of optimization methods in geostatistical studies in the past and recent years, much effort has also been made regarding the automatic fitting of the optimum model to the experimental variogram. Automatic fitting is carried out: 1) to provide a basic model to enhance fitting by eye, 2) when it is necessary that a model be fitted to a large number of experimental variograms in a short time, and 3) when it is necessary that the uncertainty related to the experimental variogram be considered in fitting a model to it (Pardo-Igúzquiza, 1999). The history goes back to formulating the problem of fitting a mel to an experimental variogram based on such criteria as least squares (David, 1977; Cressie, 1979), generalized least squares (Cressie & Hawkins, 1980, Taylor & Burugh, 1986, Genton, 1998b), weighted least squares (Starks & Fang, 1982, Cressie, 1985, Jian et al., 1996), and weighted and ordinary least squares (McBratney & Webster, 1986). In the first program developed for the automatic fitting of the variogram model, use was made of the weighted least squares criterion for the modeling of the objective function and the problem was solved using the quadratic programming method after incorporating such constraints as "monotonocity", "smoothness" and "convexity" (Shapiro & Botha, 1991). Later, use was made of the least squares criterion and the simplex optimization method in the VARFIT program (Pardo-Igúzquiza, 1999) developed in the FORTRAN environment. Genton (Gorsich and Genton 2000) made an effort in non-parametric estimation of variogram and its derivatives for the fitting of the appropriate model in the Matlab environment (version 5). Paula & Clayton (2003) enhanced the objective function of the experimental variogram considering lags intervals and the number of pairs at every point (Paula F. Larrondo et al., 2003). Desassis (Desassis & Renard, 2013) developed an algorithm, based on the numerical algorithm, to find the optimum model based on the objective function (5) for single/multi-variable cases. Many of the former methods have been developed to fit a model to a direct experimental variogram, but Emery (Emery, 2010) studied the fitting of a coregionalization model extensively and proposed 3 new algorithms for unknown sill, known sill, and Plurigaussian models. The objective function in his studies was weighted least squares and use was made of the extended Goulard-Voltz, Simulated Annealing, and Levenberg-Marqurdt methods respectively to find the optimal solutions. Effort has been made in this paper to fit the best theoretical nested variogram model to an experimental variogram. The fitted models have been selected from among spherical, exponential, gaussian, and cubic models that can be 1, 2, or 3-structured. Since most simulation algorithms and geostatistical estimations need a variogram model, this research makes possible for the geoscience researchers to select, from among the many possible cases, the most appropriate variogram model with due consideration to the number of structures. 2. Materials and methods 2.1. Problem formulation ^ For the random variable Z, the experimental variograms are calculable for different vector lags hi ; i = 1,..., n, and the optimum model should be fitted to it. An appropriate model is one ^ that its value has the least distance from those of the experimental variogram (hi) in every lag hi or, in other words, the value of the following objective function should be a minimum (David, 1977; Journel & Huijbregts, 1978; Clark, 1979): Minimize 2 ^ hi hi (2) Since the amount of the experimental variogram uncertainty is different for different lags, the weights of all the lags should not be the same (Webster & Oliver, 2001). The variogram behavior near the origin specifies the spatial continuity of the regionalized variable, and the quality of fitting a model to the primary lags is more important compared to the intermediate and end ones. Therefore, it is necessary that this point be considered in the objective function in the form of adding the value of the theoretical variogram, for the desired lag, in the denominator. The result of this correction is that the more the lag distance increases from the origin, the less its effects will be on the value of the objective function (Fig. 1). The number of existing pair points for every lag, on the basis of which the experimental variogram has been calculated, is another point to consider; the more the number of pair points, the more the reliability of the value of the experimental variogram. This point can be applied in the objective function using the coefficient of the number of existing pair points for every lag. In 1985, Cressie (Cressie, 1985) modified the above objective function in the following form to be able to consider the effects of the above 2 peculiarities: Minimize j 1 Nh ^ h j (3) ^ where Nh( j) is the number of pairs at every point in calculating the variogram, (h (j)) is the value of the experimental variogram in h( j), and (h (j); ) is the value of the theoretical variogram for the estimated parameters (). Webster& McBratney 1989 studied the effects of adding these weights and concluded at the adding improved the results. k j 1 Nh h j ^ h j (4) 0,08 0,06 0,04 0,02 0 0 10 20 30 40 50 60 70 80 90 100 Fig. 1. Trend of the reduction of the value of with respect to the lag distances for the spherical model with a sill and range of 20 and 75 respectively In fitting the optimum model, besides variogram model parameters (nugget effect, sill and range), the variogram model function £ and the number of the structures m are effective too. The rate of the effect of the variogram function increases especially in nested models where every structure can follow one of the models. In many experimental variograms, it is observed that in end lags (after reaching the value of the sill), the trends are rising, falling or sinusoidal (Fig. 3); therefore, the value of ^ ( (h (j)) (h (j); ))2 increases considerably. After range, coefficient 2 1 reaches a constant number for spherical, gaussian and cubic models (for the exponential model it tends to a constant value), and Nh( j) does not vary much in the end lag; therefore, coefficient Nh j ^ cannot compensate for the increase in ( (h (j)) (h (j); ))2 and causes the 2 h j sill of the fitted optimum model to distort. To adjust this increase, it is necessary that a "distance" function be added to the objective function. Now, since by adding h in the denominator of the objective function the weights allocated to the preliminary points will be illogically very large (compared to other points) (Fig. 2), to solve this problem, h will be substituted with h in the denominator of relation 3 and the objective function will be rewritten as follows: 0,1 0,05 0,3 0,15 Fig. 2. a) Reduction trend of h with respect to the lags distances, b) Reduction trend of h (both with respect to lags distances) Effects of the modified component (applied in the objective function) on drawing near the sill of the fitted variogram model to the real values are shown in Fig. 3. Fig. 3-a is the model fitted based on relation 3 wherein the value of the sill of the fitted variogram model is smaller than the desirable value due to the falling trend in the end lags and Fig. 3-b is related to the model fitted using the modified objective function of relation 4 wherein the value of the sill of the fitted model has increased and looks more desirable. Fig. 3. Models fitted to the experimental variogram: a) with the Cressie weighting relation and, b) with the modified relation 2.2. Decision variables Decision variables in this research are the variogram model parameters including the nugget effect, sill, and range, and are modeled in the form of the linear ((2n + 1),1) where n is the number of the structures of the variogram model. Its first array is the nugget effect and arrays 2i and 2i +1 are respectively the sill and the range of the i th structure. Another decision variable that affects the output of the algorithm is the variogram model selected for every one of the structures. The possible structures considered in the present study are spherical, gaussian, exponential, and cubic (their specifications are given in Table 1). Since the variogram model type was not definable in the form of a quantitative decision variable, the optimization problem was solved based on different arrangements of the type of model used for the variogram model structures; the better options were then selected based on the value of the objective function. TABLE 1 Specifications of the common variogram models Spherical model h h C C 3 h 2 a 1 h 2 a for h for h a, a, Exponential model Gaussian model C 1 exp 3h a 3h a C 1 exp Cubic model h h C 7 C h a 35 h 4 a 7 h 2 a 3 h 4 a for h for h a, a, 2.3. Simulated annealing The simulated annealing (SA) algorithm is a simple, effective, and meta-heuristic optimization algorithm for the solution of NP optimization problems. Similar to most meta-heuristic algorithms, it has been established by modeling and simulating one of the nature's laws or phenomena. It has been presented based on substituting physical elements in the process of physical annealing (system state, state variation energy, temperature, and freezing state) with the elements of the optimization problem (possible solution, cost, neighborhood solution, controlling parameter, and the heuristic solution) (Kirkpatrick et al., 1983). This algorithm consists of 2 basic mechanisms: 1) producing the substitute, and 2) an acceptance rule (Lee & El-Sharkawi, 2008). To solve an optimization problem, the SA algorithm first starts with a primary solution and then, in an iteration loop, moves to the neighboring solution. If the latter is better than the current one, the algorithm substitutes it (moves to it); otherwise, it will accept it as the present solution with a probability P which is found as follows: f T (5) where f is the difference between the objective function of the present solution and that of the neighboring one, and T is the temperature. Many iterations are run at every temperature and then the temperature is reduced gradually. In the preliminary steps, very high temperatures are adjusted so that worse solutions can be more probable to accept. With the gradual decrease in temperature, there is less probability for the worse solution to be accepted in the final steps; therefore, the algorithm converges on a near-optimal solution. The algorithm flow-chart is given in Fig. 4. The simulated annealing algorithm has been used to find the optimal solution of such different mining problems as optimally locating the additional drillholes (Soltani & Hezarkhani, 2009; SoltaniMohammadi & Hezarkhani, 2013), calculating the permissible charge weight for blasting operations (Soltani-Mohammadi et al., 2012), geophysics studies (Luo et al., 2013), controlling and guiding exploration and extraction operations (Debba et al., 2009; Xia et al., 2011), and simulating and estimating (Deutsch & Journel, 1992; Lee & Ellis, 1997; Peredo & Ortiz, 2011). Reasons for their high usage are: 1) they do not have to incorporate problem assumptions, and 2) they rapidly find the near optimal global solution. To run the simulated annealing algorithm, it is necessary that, first, such annealing parameters as the annealing function, the temperature updating function, and the initial temperature be specified. The annealing function is either "Boltzmann" or "fast", and the temperature updating function is selected from among exponential, logarithmic, and linear functions. The initial temperature can be defined Fig. 4. Simulated annealing algorithm both as a number or a function. 3. Case studies Three case studies have been presented to assess the capability of the proposed algorithm. 3.1. Case study 1 This case study is related to the automatic fitting of a model suitable to an omni-directional experimental variogram of an Al2O3 grade dataset with a variance equal to 21.3. This experimental variogram has been calculated in 20 lags with an initial lag of 30 m. Fig. 6 shows this experimental variogram together with the number of pair points in every lag. Running the algorithm necessitates the determination of the preliminary parameters of the annealing; to determine these parameters, different cases were studied. Fig. 5 shows the effects of the selection of the annealing function on the process of optimization. As shown, if the "fast" function is selected as the annealing function: 1) the time needed to reach the near-optimal solution will be reduced (Fig. 5c), and 2) it will be less probable for the solution to be trapped in a local optimal solution (Fig. 5(a and b)). A linear function was selected for temperature updating, and the initial annealing temperature was taken to be 300 degrees. Fig. 5. Results from running the simulated annealing algorithm: a) "fast" annealing function, b) Boltzman annealing function, and c) comparison of a) and b) regarding the best value of the objective function Fig. 6. The exponential variogram of the first dataset Fig. 7 shows the results from running the algorithm; 4 fitted, single-structured, optimum models can be seen in the first row. As shown, the spherical (7-1), cubic (7-2), and gaussian (7-3) models are fitted acceptably, but the exponential models (7-4) are not so. Four of the best 2-structured and 3-structured models fitted to the experimental variogram are given in rows 2 and 3. As shown, in Figs. (7-5) to (7-11) (except Fig. (7-6)), one structure in all the nested models is cubic. Although using the cubic model has had acceptable results in the fitting of the variogram model, it should be noted that it cannot be used in such business software brands as "Datamine", "Surpac", etc and in non-business ones such as "gslib", "Wingslib", SGeMS, etc. Cross-validation was carried out in this research based on the 12 better models (Fig. 7), and the results for Mean Error (ME), Mean Squared Error (MSE), Correlation Coefficient (CC), and Combined Error (CE) are given in Table 2. ME, MSE, and CE values are found as follows: ME 1 N 1 N Z xi i 1 N ^ Z xi (6) MSE Z xi i 1 ^ Z xi (7) Combined Error 1 abs CC MSE abs ME (8) where, for the case being studied, Z(xi) is the real value of the component (grade, thickness, etc), ^ Z(xi). is the estimated value of the component based on the fitted variogram model, and N is the number of points The least CE value is related to the single-structured gaussian model (5.1+15.6 gaussian (177)) which is shown in Fig. 7-3 as the selected model and marked with an asterisk (*) at the left corner. The worst among these 12 models is the sgle-structured exponential model. TABLE 2 CC, MSE, ME, and CE values for the first dataset n. CC MSE ME Combined Error 0.088 0.119 0.031 0.038 0.097 0.259 0.124 0.118 0.159 0.267 0.153 0.219 Fig. 7. Models proposed for the basic model of fitting by eye for the first set of data 3.2. Case study 2 Effort has been made here to automatically fit a model suitable to an Omni-directional experimental variogram of an SiO2 grade dataset with a variance equal to 8.2. This experimental variogram has been calculated in 18 lags with an initial lag length of 40 m. Fig. 8 shows the experimental variogram together with the number of pair points in every lag. Fig. 8. The experimental variogram of the second set of data Fig. 9 shows the results of running the algorithm. There are 4 fitted, single-structured, optimum models shown in the first row. As shown, the exponential model (9-1) has been fitted well, but, considering the structure of the experimental variogram, the spherical (9-2), cubic (9-3), and gaussian (9-4) models have not been fitted well. Four best 2-structured and 3-structured models fitted to the data are shown respectively in rows 2 and 3 of the figure. Similar to case study 1, in some cases, with an increase in the number of structures, the quality of the fitted model has somewhat improved, but it is negligible due to the considerable increase in the volume of calculations and the increased time required for the algorithm running and simulation. Cross-validation of this study too was carried out based on the 12 better models of Fig. 9; ME, MSE, CC, and CE results are given in Table 3. As shown, the least CE is related to the single-structured exponential model (0.14+8 exponential (189)) which is shown in Fig. 9-1 as the selected model and marked with an asterisk (*) at the left corner. TABLE 3 CC, MSE, ME, and CE values for the second set of data n. 1 CC 2 MSE 3 ME 4 Combined Error 5 0.119 0.034 0.069 0.003 0.007 0.131 0.085 0.039 0.052 0.038 -0.039 -0.021 Fig. 9. Models proposed for the basic model of fitting by eye for the second dataset 3.3. Case study 3 The experimental indicator variogram has been worked out for a thickness dataset with a variance equal to 0.22. This experimental variogram has been calculated in 20 lags with an initial lag length of 250 m, an azimuth of 135 degrees, and a dip of 0 degree. Fig. 10 shows the shape of this experimental variogram together with the number of pairs in every lag, and Fig. 11 shows the results of running the algorithm. Four fitted, single-structured optimum models are shown in the first row of Fig. 11. As shown, the cubic (11-1) and spherical (11-2) models have been fitted well, but the fitted exponential and gaussian models are not so. It is seen that the quality of the fitted model has not been improved considerably due to the increase in the number of structures. But, this small improvement will lead to the complexity of the variogram model which increases the volume of calculations and the time needed to run the krigging algorithm; therefore, it is negligible, especially when the number of grid points and data are many, or when search volume neighborhood methods are not used. Fig. 10. The experimental variogram of the third dataset Fig. 11. Models proposed for the basic model of fitting by eye for the third dataset 4. Conclusions In many experimental variograms, we observe rising, falling, or sinusoidal trends in the end lags, and this causes the variogram sill to deviate from the desirable value and go up or down; we added h in the denominator of relation 3 to modify this effect. These case studies showed that we were able, by this modification of the objective function, to minimize the effects of these trends and bring the variogram sill closer to the desirable value. As shown, increasing the number of structures from 2 to 3, does not considerably improve the quality of the fitted model, but it will lead to the complexity of the variogram model which increases the volume of calculations and the time needed to run the krigging algorithm; this will cause problems, especially when grid points and data are many, or when search volume neighborhood methods are not used. Since, in the simulated annealing, we cannot define the "equal" case when we are defining the constraints, we are not able to define the constraints related to the equality of the sill of the fitted model with the value of the data variance; we can define only the upper and lower constraints of every decision variable (co, c, a). To solve this problem, it is suggested that the future researches be focused on such other meta-heuristic algorithms as the genetic algorithm (GA).
Archives of Mining Sciences – de Gruyter
Published: Sep 1, 2016