Use of Evolutionary Computation to Improve Rock Slope Back Analysis

: Generally, in geotechnical engineering, back analyses are used to investigate uncertain parameters. Back analyses can be undertaken by considering known conditions, such as failure surfaces, displacements, and structural performances. Many geotechnical problems have irregular solution domains, with the objective function being non-convex, and may not be continuous functions. As such, a complex non-linear optimization function is typically required for most geotechnical problems to attain a better understanding of these uncertainties. Therefore, particle swarm optimization ( PSO ) and a genetic algorithm ( GA ) are utilized in this study to facilitate in back analyses mainly based on upper bound ﬁnite element limit analysis method. These approaches are part of evolutionary computation, which is appropriate for solving non-linear global optimization problems. By using these techniques with upper-bound ﬁnite element limit analysis ( UB-FELA ), two case studies showed that the results obtained are reasonable and reliable while maintaining a balance between computational time and accuracy.


Introduction
Many researchers have paid attention to assessing the stability of rock slope stability in the past few decades [1][2][3][4], but it is still a major challenge for engineers. Traditionally, design input parameters are obtained from in situ rock mechanics testing and measurement. For example, the rock mass deformation modulus can be obtained from an in situ plate loading test and the shear strength parameters from in situ block shear tests. However, due to material heterogeneity, it is challenging to accurately determine the value of strength parameters through laboratory or field tests. Back analysis can provide a very effective and rapid means of determination of these parameters.
Any failure process in rock and soil slopes can rationally be regarded as a large-scale in situ shear test performed by nature. With a proper understanding of the failure mechanism and gathering the necessary data of the failed masses, the geotechnical properties can be obtained by back analysis [5]. It can help to improve the knowledge in input parameters. In general, due to the nature of the variable, which often contains a jointed rock mass fracture, naturally occurring discontinuities and anisotropy make it a hard job for engineers to estimate the stability of a rock slope. As a result, back analysis is a very powerful method that can be implemented to determine the strength parameters of failed slopes. These back-calculated values could then be used for remedial work and the redesign of a failed slope, even new surrounding slopes with similar geotechnical conditions. Back analysis was introduced by Gioda and Maier [6] and Cividini et al. [7], for the sophistication of the observational method and establishes a crucial tool for estimating geotechnical design parameters. Back analysis has been applied in many areas in geotechnical engineering, such as: excavations [8][9][10], underground works [11,12], as well as slope stability problems [13][14][15].

Genetic Algorithm (GA)
Genetic algorithms (GA), first introduced by John Holland and his collaborators in the 1960s and 1970s, are population-based stochastic search algorithms inspired by Charles Darwin's theory of natural selection. Holland was the first to adopt the use of crossover, recombination, mutation, and selection operators in the study of artificial and adaptive systems. These operators form the basis of the genetic algorithm as a problem-solving strategy. GA begins by randomly generating a population of potential solutions. Chromosomes, which are strings of DNA, represent each member of the population, as an encoding of the problem and its search space. Members of a population are selected for recombination based on their fitness and they recombined to form a new population, for the next generation. This process is repeated until some termination criteria are met, as shown in Figure 1.
There are several control parameters in GA. The mutation probability (p m ), which is defined as the percentage of the total number of genes in the population and controls the probability with which new genes are introduced into the population. The purpose of the mutation operator is to maintain diversity within the population and suppress premature convergence. A simple way to achieve mutation would be to alter one or more genes. In GA, mutation serves the crucial role of either (a) replacing the genes lost from the population during the selection process so that they can be tried in a new context or (b) providing the genes that were not presented in the initial population. Mutation is the part of GA that is related to the "exploration" of the search space.
The crossover operator plays an important role in the overall operation of GA. The performance of GA to a great extend depends on the performance of the crossover operator used. It is analogous to reproduction and biological crossover. In this, more than one parent is selected and one or more offspring are produced using the genetic material of the parents. The crossover probability (P C ) is defined as the probability of the number of offspring produced in each generation of the population size. Crossover is usually applied in a GA with a high probability.
Genetic algorithms are able to address complicated problems with many variables and a large number of possible outcomes by simulating the evolutionary process of "survival of the fittest" to reach a defined goal. They operate by producing many random answers to a problem, eliminating the worst and cross-pollinating better answers. This process of repeated elimination and regeneration gradually improves the quality of the answers to an optimal or near-optimal condition. Appl. Sci. 2019, 9,

Particle Swarm Optimization (PSO)
Particle swarm optimization (PSO), inspired by the flocking and schooling patterns of birds and fish, was first invented by Russel Eberhart and James Kennedy in 1995 [25]. A group of variables has its values altered to get close to the member whose value is closest to the targeted solution several iterations at any given time. Depict a flock of birds circling over an area where they can sense a hidden food source. The bird that is closest to the food chirps the loudest and the other birds swing around in his direction. If one of the other birds comes closer to the target than the previous one, it chirps louder and the others sway over to him. This pattern continues until one of the birds find the food. The whole concept of PSO is based on this.
First, the algorithm is initialized by generating a random number of particles and a fitness function is used to evaluate the fitness of each particle. Each particle has a velocity, which is used to direct its flying pattern. The flying of each particle is altered by its own flying experience and those of its companions' flying experience this scenario is a good representation of a possible solution to a problem. The best solution obtained by a particle so far is known as its personal best (Pbest) and the best solution obtained by any of the particles in the swarm so far is known as the global best ( ).

Particle Swarm Optimization (PSO)
Particle swarm optimization (PSO), inspired by the flocking and schooling patterns of birds and fish, was first invented by Russel Eberhart and James Kennedy in 1995 [24]. A group of variables has its values altered to get close to the member whose value is closest to the targeted solution several iterations at any given time. Depict a flock of birds circling over an area where they can sense a hidden food source. The bird that is closest to the food chirps the loudest and the other birds swing around in his direction. If one of the other birds comes closer to the target than the previous one, it chirps louder and the others sway over to him. This pattern continues until one of the birds find the food. The whole concept of PSO is based on this.
First, the algorithm is initialized by generating a random number of particles and a fitness function is used to evaluate the fitness of each particle. Each particle has a velocity, which is used to direct its flying pattern. The flying of each particle is altered by its own flying experience and those of its companions' flying experience this scenario is a good representation of a possible solution to a problem. The best solution obtained by a particle so far is known as its personal best (Pbest) and the best solution obtained by any of the particles in the swarm so far is known as the global best (Gbest). The velocity and position of each particle are updated using its Pbest and Gbest according to the equations below. Velocity: where, vk(t): Velocity of the particle k at iteration t c 1 = c 2 : Constants w: Inertia weight. r 1 = r 2 : Random numbers ranging for [0, 1] Pbest: Personal best Gbest: Global best xk(t): Current position of particle k at iteration t Position: Equation (1) is used to determine the new velocity of particle k. The inertia weight (w) determines the extent to which the particle remains along its original course unaffected by the pulls of P best and G best [25]. c 1 and c 2 imitate the weights whether returning to the location of the personal best or exploring the location toward the global best and the values can be determine through empirical study. The purpose of r 1 and r 2 is to mimic the unpredictable behavior of natural swarms. Particle k flies towards a new position according to Equation (2). A fitness function, which is usually the objective function under consideration, is used to measure the performance of each particle.

Optimization Strategy for Back Analysis
Many geotechnical and other engineering problems can be formulated as optimization problems, for example stability-type problems. An optimization algorithm is a useful and powerful technique in both deterministic and probabilistic slope stability analysis.
Generally, optimization algorithms can be divided into two basic classes: deterministic and heuristic/metaheuristic algorithms. The deterministic methods are based on a gradient function. However, the acquisition of gradient information can be very difficult or, in some cases, altogether impossible to obtain. The metaheuristic method is not restricted in this manner. The metaheuristic methods include algorithms such as genetic algorithms, particle swarm optimization, simulated annealing, Tabu search, and ant colony optimization. The metaheuristic methods are derivative free methods and can be applied to any optimization problem irrespective of the linearity or non-linearity of the objective function to be minimized and constraints.
The application of optimization in slope stability problems was first investigated by Baker in 1980 [26]. In this approach, Baker combined an optimization search using the dynamic programming technique with Spencer's method [27] to locate the critical slip surface and calculate the associated factor of safety. The method proposed by Baker can be applied to slopes of any geometry, layering, pore pressure and external load distribution and could find a more critical slip surface with a lower factor of safety compared to preceding methods.
Other studies [28,29] proposed the genetic algorithm methodology for determining the critical slip surface for circular and multiple-wedge stability analysis, respectively. The results of these two studies demonstrated the effectiveness of the genetic algorithm approach. Figure 2 shows the optimization strategy for back analysis. The numerical model simulates the construction process and the geometry of geotechnical structures. Then, the optimization algorithm performs an iterative process computing solution for the geotechnical parameters to minimize the objective function. The objective function also known as the error function assesses the discrepancy between the computed and the monitored or measured values.

Objective Function
The solution of the objective function is based on some optimization techniques for determining a set of input parameters that make the value of the objective function minimum [31]. In this study, the objective function involves using the summed squared errors between the measured mean value unit weight under which the slope collapse ( ) and the determined upper-bound finite element limit analysis (UB-FELA) unit weight ( ). The essence of the objective function is to minimize the difference between ( ) and ( ) by successively fitting the simulation result to the measured mean value unit weight under which the slope failed so as to determine the uncertainty parameters of the rock slope that are of interest to us.
where is the objective function to be minimized.

Objective Function
The solution of the objective function is based on some optimization techniques for determining a set of input parameters that make the value of the objective function minimum [30]. In this study, the objective function involves using the summed squared errors between the measured mean value unit weight under which the slope collapse (γ m ) and the determined upper-bound finite element limit analysis (UB-FELA) unit weight (γ F ). The essence of the objective function is to minimize the difference between (γ m ) and (γ F ) by successively fitting the simulation result to the measured mean value unit weight under which the slope failed so as to determine the uncertainty parameters of the rock slope that are of interest to us.
where δ is the objective function to be minimized.

Overview of Slope Failure at a Barite Open Pit Mine
The Baskoyak mine in central Turkey is an open mine pit operated for the extraction of barite. A comprehensive slope stability project was conducted to determine the engineering properties of the rock mass, and to assess the failure mechanism and the alternatives for improving the overall stability between 1987 and 1988, and the investigation was published by Ulusay and Aksoy [31]. Based on the scanline surveys, Ulusay and Aksoy reported that the schist should be regarded as comprising of two rock mass types. The first type comprises a schist rock mass, heavily broken by closely spaced discontinuities and the second type is a weathered schist in different degrees both in the hanging wall and footwall. The unit weight of the schists ranges between 17.2 KN/m 3  The rock mass rating (RMR) value of the rock mass was calculated using the 1989 version of RMR classification [32] using parameters given by Sonmez et al. [33]. As indicated by Sonmez and Ulusay [34], no sign of groundwater was encountered. Thus, the pit slopes were treated as dry for stability assessments.

Numerical Model of the Rock Slope
The upper bound finite element limit analysis method was used in conjunction with the latest version of the Hoek Brown failure criterion [35] to determine the unit weight of the failed rocked slope. Figure 3 shows the model geometry of the slope cross-section considered for analysis. Table 1 shows the input parameters employed in the analysis.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 21 scanline surveys, Ulusay and Aksoy reported that the schist should be regarded as comprising of two rock mass types. The first type comprises a schist rock mass, heavily broken by closely spaced discontinuities and the second type is a weathered schist in different degrees both in the hanging wall and footwall. The unit weight of the schists ranges between 17.2 KN/m 3  The rock mass rating (RMR) value of the rock mass was calculated using the 1989 version of RMR classification [33] using parameters given by Sonmez et al. [34]. As indicated by Sonmez and Ulusay [35], no sign of groundwater was encountered. Thus, the pit slopes were treated as dry for stability assessments.

Numerical Model of the Rock Slope
The upper bound finite element limit analysis method was used in conjunction with the latest version of the Hoek Brown failure criterion [36] to determine the unit weight of the failed rocked slope. Figure 3 shows the model geometry of the slope cross-section considered for analysis. Table 1 shows the input parameters employed in the analysis.

GA Parameter Settings
There are many parameter and operator settings in genetic algorithm that can affect the performance of the algorithm. Such settings can be implemented differently in various problems. Table 2 shows the parameter settings used in this study based on the recommendations given by DeJong et al. [37].  There are many parameter and operator settings in genetic algorithm that can affect the performance of the algorithm. Such settings can be implemented differently in various problems. Table 2 shows the parameter settings used in this study based on the recommendations given by DeJong et al. [36]. The choice of PSO parameters can have a significant influence on optimization performance. As a result, selecting PSO parameters that yield good performance has been the subject of much research. In this study, the parameter settings of PSO were done in accordance with the recommendations given by Shi et al. [25]. Table 3 shows the parameter settings of PSO used throughout this study. Herein the disturbance factor (D) of the failed rock slope is back analyzed. The uncertainty parameter D is back calculated considering a uniform distribution and a range of 0-1. Back calculation is accomplished by using the optimization algorithms to minimize the objective function (Equation (3)). Table 4 shows the back calculated D values obtained from PSO and GA; whereas Figures 4 and 5 show the convergence characteristic of the objective function from PSO and GA, respectively.     Both PSO and GA yield the same back calculated D value of 0.68 as shown in Table 3, with a resulting predicted unit weight of 3 The convergence characteristics of the objective function in Figures 5 and 6 shows that the resulting predicted unit weight perfectly matches the measured mean value unit weight under which the slope collapse

Back Calculation of Two Uncertainty Slope Parameters
This involves considering the back calculation of two uncertainty parameters D and geological strength index (GSI) of the rock slope. The parameters D and GSI are back calculated considering uniform distributions with ranges of 0-1 and 11-21 respectively. Tables 5 and 6 respectively show the back calculated D and GSI values obtained from PSO and GA. Both PSO and GA yield the same back calculated D value of 0.68 as shown in Table 3, with a resulting predicted unit weight of γ F = 22.2 KN/m 3 The convergence characteristics of the objective function in Figures 5 and 6 shows that the resulting predicted unit weight perfectly matches the measured mean value unit weight under which the slope collapse γ m = 22.2 KN/m 3 .  Figures 6 and 7 shows that the resulting perfectly matches the .     Tables 5 and 6 closely match the back calculated D value from the single parameter back calculation as shown in Table 4.   In the single parameter back calculation, the same back calculated D value was obtained for both PSO and GA. In addition, the double parameters back calculation yield very similar back-calculated values obtained from both PSO and GA. However, in the three parameters' back calculation, the back-calculated values of the parameters D and m i from PSO and GA are quite different as shown in Tables 7 and 8. In addition, the back-calculated values of D and GSI in this case are quite different from those obtained in Tables 4 and 5, for the two uncertainty parameters back calculation.

Back Calculation of Four Uncertainty Slope Parameters
Herein, the uncertainty parameters D, GSI, m i and σ ci of the rock slope are considered for back calculation. The uncertainty parameters D, GSI, m i and σ ci were all back-calculated considering uniform distributions with ranges of 0-1 for D, 11-21 for GSI, 4-10 for m i and 2-8 MPa for σ ci . Tables 9 and 10, respectively, show the back calculated D, GSI, m i and σ ci values obtained from PSO and GA. Back calculation is accomplished by using the optimization algorithms to minimize the objective function (Equation (3)).
Tables 9 and 10 show that the back-calculated values of all the parameters from PSO and GA are quite different, as experienced in the three parameters' back calculation. This is in contrast with the single parameter back calculation, which yields the same value of D, and double parameter back calculation, which yield very similar back-calculated values of D and GSI from PSO and GA. This shows that as the number of uncertainty parameters increase, different methods of back calculation will yield different back calculated values. The convergence characteristics of the objective function in Figures 6 and 7 shows that the resulting γ F perfectly matches the γ m .

Comparison of Failure Mechanisms
To ascertain the validity of the back calculated parameters, failure mechanisms of the backcalculated parameters from one, two, three and four uncertainty parameters back calculations for both PSO and GA are compared with the original failure mechanism of the slope. The gray colors shown within the slope cross-section in Figure 8 show the nature of the failure mechanisms of the slope. Figure 8 shows satisfactory similarities in the failure mechanisms between back calculated values from PSO and GA compared with the original slope failure for one and two uncertainty slope parameters. The same phenomenon has been observed for back calculation of three and four uncertainty slope parameters. Therefore, it can be deduced that for a reasonable range, values obtained from back calculation do not influence the failure mechanism significantly.

Comparison of Failure Mechanisms
To ascertain the validity of the back calculated parameters, failure mechanisms of the backcalculated parameters from one, two, three and four uncertainty parameters back calculations for both PSO and GA are compared with the original failure mechanism of the slope. The gray colors shown within the slope cross-section in Figure 8 show the nature of the failure mechanisms of the slope.  Figure 8 shows satisfactory similarities in the failure mechanisms between back calculated values from PSO and GA compared with the original slope failure for one and two uncertainty slope parameters. The same phenomenon has been observed for back calculation of three and four uncertainty slope parameters. Therefore, it can be deduced that for a reasonable range, values obtained from back calculation do not influence the failure mechanism significantly.

Statistical Study on the Distribution and Range of Disturbance Factor (D)
As shown in the failure surfaces from one, two, three and four uncertainty parameters, the back calculations failure mechanism is almost unchanged. However, it is important to obtain accurate back calculation results in practice. A study was conducted to determine how the mode of distribution and the range chosen for back calculating D influence the other back calculated parameters. A statistical study was conducted considering both normal and uniform distributions to evaluate how the distribution and range of D influence the coefficients of variation (COVs) of the other back-calculated parameters.
The range and distribution for , and are based on the recommendations given by Hoek. From the experience of Hoek , and were normally distributed with COVs of 0.25, 0.125 and 0.1 respectively. The original values given by Sonmez et al. [34] were chosen as the average values in this study as shown in Table 1. Based on the above recommendations the standard deviation and range of the parameters , and are determined for back analysis as shown in Table  11.
Original failure mode Back calculating D from PSO and GA Back calculating D and GSI from PSO Back calculating D and GSI from GA

Statistical Study on the Distribution and Range of Disturbance Factor (D)
As shown in the failure surfaces from one, two, three and four uncertainty parameters, the back calculations failure mechanism is almost unchanged. However, it is important to obtain accurate back calculation results in practice. A study was conducted to determine how the mode of distribution and the range chosen for back calculating D influence the other back calculated parameters. A statistical study was conducted considering both normal and uniform distributions to evaluate how the distribution and range of D influence the coefficients of variation (COVs) of the other back-calculated parameters.
The range and distribution for σ ci , m i and GSI are based on the recommendations given by Hoek. From the experience of Hoek σ ci , m i and GSI were normally distributed with COVs of 0.25, 0.125 and 0.1 respectively. The original values given by Sonmez et al. [33] were chosen as the average values in this study as shown in Table 1. Based on the above recommendations the standard deviation and range of the parameters σ ci , m i and GSI are determined for back analysis as shown in Table 11. Since there is no given recommendation on the COV and mode of distribution for D, this study uses different approaches to find a scenario, which could be employed for back calculating D.
Scenario A: Considering a uniform distribution and a range of 0-1 for back analyzing the disturbance factor (D).
Scenario B: Considering a normal distribution and a range of 0-1 for back analyzing the disturbance factor (D).
Scenario C: Considering a normal distribution and a range of 0.4-1 for back analyzing disturbance factor. This aims to evaluate the impact of reducing the distributed range of D on the other back-calculated parameters. In fact, the range of D can be smaller if the slope surface condition is judged by experienced engineers. In addition, the study by Li et al. [37] indicated D = 0.7 would be reasonable for this case. Therefore, D = 0.7 is assumed as the mean value herein.
For each of the three different scenarios, all four parameters (D, GSI, m i and σ ci ) are back calculated 30 times so as to determine the COVs of all four parameters for all three different scenarios. Table 12 showed the mean values obtained from PSO and GA. While Table 13 showed the comparison of COVs obtained from PSO and GA from the statistical study. The mean value for each parameter is determined for both PSO and GA for each of the three scenarios taken into consideration for the statistical study, as shown in Table 12. The mean values of the back calculated parameters are close to the original strength parameters in each of the scenarios taken into consideration for both PSO and GA. For the single and double parameter back calculation, the values of the back-calculated parameters closely match those of the mean values from the statistical study. However, as the number of uncertainty parameters increase to three and four parameters the back-calculated parameters differ from the mean values obtain from the statistical study significantly.  Table 13 shows the comparison of COV values for PSO and GA obtained from the statistical study. The COVs obtained from GA are smaller when compared to PSO for all four parameters in all three scenarios. In addition, it can be observed that COVs of GSI and m i are not significantly influenced by the mode of distribution and range of D considered for back calculation. The results also indicate that if distributed range of D can be reduced, the COV for back-calculated σ ci decreases noticeably. Overall using a normal distribution gives better results compared to uniform one in terms of accuracy.

Introduction
Strong earthquakes can result in large amounts of landslide thus resulting in human casualties and property damage. In recent years, there have been many reports of serious damages caused by earthquake-induced landslides [38]. The 2008 Wenchuan earthquake shocked the Sichuan province in china and caused as many as 60,104 landslides [39]. The Daguangbao landslide, with an approximate influenced area of 7.3-10 million m 2 and a volume of 750-850 million m 3 , was the biggest landslide induced by the 2008 Wenchuan earthquake. Not only was it the largest know landslide in China, but also one of the largest known landslides all over the world. Several studies conducted on the Daguangbao landslide [40] focused on the characteristics, failure mechanism, and geological properties immediately after sliding. In addition, several post-earthquake field investigations were also carried out [41], to further study the characteristics and failure mechanism.

Physical and Mechanical Properties of the Rock Slope
According to Zhang et al. [42], experience equations based on that of Hoek-Brown [43] were used to determine the strength parameters of the slope. Experience values of the input index for the Hoek-Brown criterion are listed in Table 13. Zhang et al. [42] indicated that there is some degree of uncertainties associated with the estimation of material properties such as GSI, m i , and D. The unit weight of the slope is determined using the upper-bound finite element limit analysis (UB-FELA) approach using the parameters given in Table 14 as input parameters. As indicated by Marinos et al. [44], rock mass disturbance (D) should be considered because of its significant influence on the evaluation of rock slope stability. The damage zone thickness (T) refers to the thickness of the disturbed zone of the slope. Figure 9 shows a rock slope divided into two regions to account for the disturbed zone.
in china and caused as many as 60,104 landslides [41]. The Daguangbao landslide, with an approximate influenced area of 7.3-10 million m 2 and a volume of 750-850 million m 3 , was the biggest landslide induced by the 2008 Wenchuan earthquake. Not only was it the largest know landslide in China, but also one of the largest known landslides all over the world. Several studies conducted on the Daguangbao landslide [42] focused on the characteristics, failure mechanism, and geological properties immediately after sliding. In addition, several post-earthquake field investigations were also carried out [43], to further study the characteristics and failure mechanism.

Physical and Mechanical Properties of the Rock Slope
According to Zhang et al. [44], experience equations based on that of  were used to determine the strength parameters of the slope. Experience values of the input index for the Hoek-Brown criterion are listed in Table 13. Zhang et al. [44] indicated that there is some degree of uncertainties associated with the estimation of material properties such as GSI, mi, and D. The unit weight of the slope is determined using the upper-bound finite element limit analysis (UB-FELA) approach using the parameters given in Table 14 as input parameters. , rock mass disturbance (D) should be considered because of its significant influence on the evaluation of rock slope stability. The damage zone thickness (T) refers to the thickness of the disturbed zone of the slope. Figure 9 shows a rock slope divided into two regions to account for the disturbed zone.

Geometry of the Upper-Bound Finite Element Limit Analysis (UB-FELA) Model
The main sliding direction of the Daguangbao landslide is chosen as the slope profile for the UB-FELA model as shown in Figure 10. For the simulation, the whole slope is split into three parts: the base block, the upper sliding mass, and the lower sliding mass. The splitting is undertaken based on the failure surface and topography of the slope. The UB-FELA simulation is conducted considering different thickness to height ratios ( / ) for the slope. / = 2.5 and / = are the two thickness to height ratios considered in the UB-FELA simulation subjectively.

Geometry of the Upper-Bound Finite Element Limit Analysis (UB-FELA) Model
The main sliding direction of the Daguangbao landslide is chosen as the slope profile for the UB-FELA model as shown in Figure 10. For the simulation, the whole slope is split into three parts: the base block, the upper sliding mass, and the lower sliding mass. The splitting is undertaken based on the failure surface and topography of the slope. The UB-FELA simulation is conducted considering different thickness to height ratios (T/H) for the slope. T/H = 2.5 and T/H = in f inity are the two thickness to height ratios considered in the UB-FELA simulation subjectively.   The upper bound finite element limit analysis is used in conjunction with the Hoek-Brown failure criterion.

Back Analyzing the Horizontal Seismic Coefficient (K h ) Slope Parameter
In this study, using the optimization algorithms (PSO and GA) and Equation (3) as the fitness function, the unit weight (γ F ) obtained from the UB-FELA method is fitted successively to match the targeted unit (γ m ) by back analyzing for the horizontal seismic coefficient (k h ) parameter of the rock slope.
Herein, the back calculated horizontal seismic coefficient is the critical k h value of the slope. As shown in Table 15 for T/H = 2.5 the critical k h = 0.120 and for T/H = infinity, the critical k h = 0.045. For the same thickness to height ratio, the values of the back-calculated k h parameters are the same for both PSO and GA. The results also show the importance of considering thickness to height ratio in rock slope stability analysis. The critical k h value is much smaller for T/H = infinity compared to T/H = 2.5. The uncertainty parameter σ ci of the upper and the lower sliding masses of the slope as shown in Figure 10 as well as the horizontal seismic coefficient (k h ) of the slope are considered for back calculation. Back calculation is accomplished by using the optimization algorithms (PSO and GA) to minimize the objective function (Equation (3)). Tables 16 and 17  Herein, the uncertainty parameters σ ci and GSI of the lower and the upper sliding masses as shown in Figure 10 as well as the k h parameter of the entire slope are considered for back calculation. GSI is considered to be same for both the upper and the lower sliding masses. In his study, Sonmez et al. [33] gave the same GSI value for the upper and the lower sliding masses. Marinos et al. [44] indicated that the GSI system works well for engineering geologists since it is consistent with their experience in describing rocks and rock masses during logging and mapping. Therefore, GSI for material 2 and material 3 are considered to be the same during back calculation. Tables 18 and 19 respectively showed the back-calculated parameters of the slope obtained from GA and PSO. The back calculated σ ci and GSI values of the lower and upper sliding masses of the rock slope obtained from the different optimization methods are different. Also the value of k h obtained from the different optimization methods are different for the same thickness to height ratios. Table 18 shows the same k h values considering T/H = infinity and T/H = 2.5 which is in contrast with phenomenon observed from the other tables where k h values are always larger for T/H = 2.5. Considering parameter distributions in this case can help us to capture more reasonable strength parameters. The uncertainty parameters σ ci , GSI and m i of the lower and the upper sliding masses as shown in Figure 10, as well as the k h parameter of the entire slope are considered for back calculation. Tables 20  and 21    The light green colors shown within the slope cross-sections in Figures 11 and 12 show the nature of the failure mechanisms. The results show satisfactory similarities in the failure mechanisms for parameters obtained from back calculation and the original slope failure mode. This indicates that with reasonable ranges for back calculation the failure mode is not significantly affected by the values of the parameters obtained from back calculation. In addition, T/H = 2.5 give a better representation The light green colors shown within the slope cross-sections in Figures 11 and 12 show the nature of the failure mechanisms. The results show satisfactory similarities in the failure mechanisms for parameters obtained from back calculation and the original slope failure mode. This indicates that with reasonable ranges for back calculation the failure mode is not significantly affected by the values of the parameters obtained from back calculation. In addition, T/H = 2.5 give a better representation of the failure mode when compared to the original failure mode as indicated by the dark line drawn on the original failure mode shown in Figure 12 when compared to that shown in Figure 11.

Statistical Study
Herein, a statistical study is conducted using both optimization algorithms to back calculate the uncertainty parameters σ ci , GSI, m i and k h while considering different thickness to height ratios. The aforementioned parameters are back calculated 30 times considering T/H = infinity and T/H = 2.5. The range for σ ci , m i and GSI are based on the recommendations given by Hoek. From the experience of Hoek σ ci , m i and GSI were normally distributed with COVs of 0.25, 0.125 and 0.1, respectively. The material properties shown in Table 14 were chosen as the mean values for back calculation. The range chosen for back calculating each uncertainty parameter is shown in Table 22. The study considers a uniform distribution for back calculating the uncertainty parameters. The parameter k h is the same throughout the slope. The uncertainty parameter GSI is the same for both the upper and the lower sliding mass. However, σ ci and m i for the upper and the lower sliding masses are different. As a result, σ ci and m i for the upper sliding mass are denoted as σ ci (3) and m i (3) and σ ci (2) and m i (2) for the lower sliding mass. Tables 23 and 24 show the comparisons of the mean values when T/H = 2.5 and T/H = infinity respectively for both PSO and GA. As expected, k h for T/H = infinity is smaller than that of T/H = 2.5. However, for GA results, k h for T/H = infinity is very similar to that for T/H = 2.5. In Table 24, it can be seen that σ ci for Material 2 is much larger than the given average value. This revealed that σ ci has significant effects on rock slope stability. While Tables 25 and 26 show the comparisons of COVs for T/H = 2.5 and T/H = infinity respectively. Results from both tables show that the COVs for k h and GSI are quite close for the different optimization algorithms as well as for different T/H ratios.

Conclusions
Procedures for rock slope back calculation considering material uncertainties were presented using as an example a slope failure in an open mine pit in Turkey as case study 1 and a slope failure induced by an earthquake in China as case study 2. The disturbance factor (D), geological strength index (GSI), jointed rock mass (m i ) and unconfined compressive strength (σ ci ) are determined to be key parameters in case study 1. Case study 2 involves seismic activity as a result of the parameters considered for back calculation being the horizontal seismic coefficient (k h ), GSI, jointed rock mass (m i ) and the unconfined compressive strength (σ ci ) which are determined to be key parameters. It was observed that as the number of uncertainty parameters increase, the results of the back-calculated values obtained from PSO and GA become very different.
In this study, using upper-bound finite element limit analysis in conjunction with PSO or GA, back analyses were successfully developed and performed. The two procedures for back calculation yield results that are identical for the single parameter back calculation.
Overall, a normal distribution gives better results compared to a uniform distribution in term of accuracy based on the COV values obtained from the statistical studies.
Both approaches of back calculation are easy to implement and require little computational effort and time. In general, for the same population size and number of iterations, the computational time for GA is 3-5 times less than that of PSO. For case study 1, GA gave smaller COV values compared to PSO but this phenomenon is not repeated in case study 2. As a result, such a phenomenon might be case-specific.
Failure mechanisms are not changed with different back-calculated strength parameters. However, parameter distributions can still help us to capture more reasonable strength parameters. This can enhance our knowledge of input parameters for the slope stability assessment or further remediation methods' selection.
By comparing the execution time in this study, PSO is about 3 times slower than GA, and it can be concluded that genetic algorithms are superior to particle swarm optimization for the back analysis problem, at least concerning rock slopes.