Back Analysis of Geomechanical Parameters Using Hybrid Algorithm Based on Difference Evolution and Extreme Learning Machine

Since the geological bodies where tunnels are located have uncertain and complex characteristics, the inverse problem plays an important role in geotechnical engineering. In order to improve the accuracy and speed of surrounding rock identification, the back analysis objective function with usage of the displacement and stress monitoring data has been constructed, with a hybrid algorithmproposed. An extreme learningmachine (ELM) is employedwith optimal architecture trained by the difference evolution (DE) arithmetic. First, the three-dimensional numerical simulation is used in the creation of training and testing samples for ELM model construction. Second, the nonlinear relationship between rock parameters and displacement is constructed by numerical simulation. Finally, the geophysics parameters are obtained byDEoptimization arithmetic taking into consideration themonitoring data including both displacement and pressure. This method had been applied in the Fusong highway tunnel in Fusong City of China’s Jilin Province, with a good effect obtained. It takes full advantage of DE and ELM and has both calculation speed and precision in the back analysis.


Introduction
In the inverse problem, the aim is to reconstruct the model or identify the parameters from a set of measurements.Since the geological body has both uncertainty and complexity, the inverse problem plays an important role in geotechnical engineering.From the degree of mathematics, the back analysis includes exact analytical method, half analytical method, and numerical method [1].Sakurai et al. proposed back analysis thinking for rock parameters data [2,3].In addition, Mashimo proposed the concept of "redesign" by summarizing the updated tunnel engineering techniques in Japan [4].Moreover, Asadollahpour et al. studied the back analysis of closure parameters of the Panet equation and Burger's model for Babolak water tunnel conveyance [5].More and more scholars paid attention to tunnel back analysis research in recent years [6][7][8].Gioda et al. studied and compared the back analysis procedures for the interpretation of field measurements in geomechanics [9][10][11].Baroth and Malecot presented the probabilistic analysis of inversely analyzing an excavation problem [12].
Despite the mathematical elegance of the exact nonlinear inversion schemes, they are limited in applicability, for the exact inversion techniques are usually applicable only in idealistic situations that may not hold in practice.For example, even a very normally true Horseshoe-Shaped tunnel with supporting cannot be expressed by exact analytical equations.Because of the limitation in monitoring techniques and the difficulty in analysis, the measured displacement is only considered as input information in the conventional back analysis, which is called "displacement back analysis" in geotechnical engineering.Along with the development of the monitoring means, there is a variety of monitoring information that can be collected by long distance automation transmission systems with the sensors of earth pressure cell, steel bar meter, strain gauge, and so forth.So the new back analysis method with usage of more information needs to be developed [13].
There is the nonlinear and implicit mapping relation between the parameters and the observation indices of geotechnical engineering, so how to express it becomes the study focus.The direct and effect model is the numerical expression such as Finite Element Method, which is based on the mechanics theory; however it is limited because the numerical model consumes too much time in back analysis process [14].Polynomial expression is easy to use for the expression, but its reasonable structure is difficult to construct; therefore the machine learning methods such as Artificial Neural Network (ANN) and support vector machine (SVM) are increasingly used in the expression of geotechnical model [15].Wu et al. respectively applied SVM in surrounding rock displacement forecast and analyzed the kernel function and parameters affecting forecast error [16,17].
Back analysis is an optimization process in essence, with the steepest descent method, the conjugate gradient method, inexact line search, and trust region algorithm [18] used generally in conventional back analysis.Because geotechnical engineering is complex and nonlinear problem cannot be easily expressed by explicit formula, the derivation cannot be realized by the above-mentioned conventional optimization methods.Furthermore, they are easily limited in local optimum solution.Therefore the soft computing methods are used in the engineering back analysis.Cui and Sheng applied genetic algorithms (GA) in the probabilistic finite element analysis of geotechnical problems [19].Karaboga et al. studied artificial bee colony (ABC) algorithm in the engineering, respectively [20,21].Jiang et al. used the particle swarm optimization (PSO) algorithm in tunnel engineering back analysis [22,23].Although the use of the tunnel analysis technologies has yielded many beneficial researching results, there are also some problems in the following: (1) Back analysis methods adopted only displacements as input data, and the computing models are often simplified to 2D patterns that are incapable of reflecting the three-dimensional space effect of a tunnel; thus, the method taking usage of three-dimensional numerical model and more information is necessary.
(2) The parameters back analysis based on a threedimensional numerical model cost too much time, and the ANN topological structure is too complex, with its completion depending on artificial experience.And the conventional optimization method is easily limited in a local optimum solution.
(3) Because the model complexity and measuring indices are inconsistent, "displacement back analysis" cannot continue to be used.Genetic algorithms have complex operation processes, and particle swarm optimization converges irregularly for having no strict convergence theory background.
Therefore, the paper considered the new excellent soft computing algorithms for difference evolution (DE) arithmetic and extreme learning machine (ELM) and proposed a new hybrid algorithm combining ELM with DE, adapting to the geomechanical parameters back analysis of multiinformation (including geophysical prospecting, displacement, and pressure data) and the three-dimensional numerical model.It is also applied in a highway tunnel in Fusong city in China's Jilin Province.

The Tunnel Back Analysis of Hybrid Algorithm Based on DE-ELM
2.1.The Inverse Problem of Geotechnical Engineering.The rock parameters should be back analyzed by field monitoring, which is opposite to the forward analysis (Figure 1).Back analysis includes the analytic and numerical methods.The numerical method theory is as follows: Select a group of parameters and calculate the surrounding rock deformation or stress; if there is a difference between calculation results and monitoring data, the input parameters will be adjusted and recalculated till the difference between the calculation and monitoring is small enough; the corresponding parameters are the back analyzed results.
In the calculation process (Figure 1), constructing the correct nonlinear relation between the surrounding rock parameters and surrounding rock displacement (or stress) remains the key issue, so we use the DE-ELM model to express the relation.It constructs the tunnel back analysis mathematic model in the following: where  (z 1 , z 2 , . . ., z r ) is the  rock mechanics parameters (or stress) to be back analyzed.Their scopes are determined by the above rock mass classification results, ( ) is the implicit function related to displacements and rock mechanics parameters,  0  is monitoring displacements (or pressure, stress) of tunnel-surrounding rock,   is calculation displacements of tunnel-surrounding rock, and  is the number of monitoring points.

The Difference Evolution
Algorithm.The differential evolution (DE) algorithm, proposed by Rainer Storn and Ken.Price, is a new type of directly global optimization algorithm that surfaced in 1995 [24].DE algorithm had NP parameter vectors from generation G Newly generated parameter vector been successfully used in optimizing ANN model and tunnel parameter identification [25,26].The specific process is as follows.
(1) Generate Initial Population.DE's most versatile implementation maintains a pair of vector populations, both of which contain   -dimensional vectors of real-valued parameters.The solution vector can be expressed as  , = ( 1 ,  2 , . . .,   ),  = 1, 2, . . .,   .  is the number of solution vectors in the  generation;  is the generation of evolution population and  is the location of an individual in the population.Consider where rand is a random number within [0, 1] and  ,1 is the solution vectors of the first generation.  is the number of solution vectors, and    and    are the upper and lower bounds of the variables, respectively.
(2) Mutation.DE generates new parameter vectors by adding the weighted difference between two population vectors to a third vector in an operation called mutation.As for the target vectors in the  generation, each one contains the component .The mutated vector is where  1 ,  2 , and  3 are a random integer from [1, 2, . . .,   ], which is different from the others.V ,+1 is the mutated vector, and  is one of the main parameters used to control the amplification of the differential variation and is a mutagenic factor whose value is taken between 0.5 and 1.The generating process of the mutated vector V ,+1 in the two-dimension solution space is shown in Figure 2.
(3) Crossover.The mutated vector's parameters are mixed with the parameters of another predetermined target vector according to certain rules in order to generate a trial vector, which is then called "crossover." The target vector  , and the mutated vector V , cross based on the following rules are to generate a trial vector  ,+1 : where where rand () ∈ [0, 1] is the component  of the vector corresponding to random number, CR ∈ [0, 1] is the crossover factor, and () is an integer randomly selected in 1, 2, . . ., .
(4) Selection.The greedy search method used to select the new solution vector  ,+1 is DE algorithm offers several significant advantages over some widely used traditional iteration algorithms such as trust region algorithms (TRA) and genetic algorithm (GA).TRA needs to solve the subproblem based on descent gradient of the objective function at each iteration step, so it is not suitable for implicit objective functions without explicit expression.GA applies binary coding and has single crossover operation with the mutation type being scant.It also has many parameters affecting the optimization ability.However, DE encodes all parameters as floating-point numbers, regardless of their type.Even integer and discrete variables are encoded as real values to add diversity to their difference distributions.DE mutation operation adopts the different strategy, that is, the adoption of the difference vector of the individuals in the population to disturb the individual, for the individual variation to be realized.The above mutation type of DE takes usage of the population distribution characteristics, improving the algorithm searching ability.And DE has only two control parameters, making it more easy to use.

The Back Analysis Method Based on DE-ELM Hybrid
Algorithm.Extreme learning machine (ELM) is the single hidden layer feed forward neural networks (SFLNs) that was proposed in 2005 after BP networks and the support vector machine [27].Its advantage is that it has simple structure, strong generalization ability, and quick learning speed, allowing it to avoid the problems of local minimum and too many iterations.For  different learning samples (  ,   ) ∈   ×   ( = 1, 2, . . ., ), the th sample output value of SLFN with  hidden layer nodes and a hidden layer activation function of () is as follows: where It can be simply expressed as where where   ,   , and   have the same meaning as formula (7),  is the output array of the hidden layer of the neural network, (  ) is the th row vector, and the th column of  array is the output corresponding to the th hidden layer node while the input variables are  1 ,  2 , . . .,   .Aiming at certain training samples, the standard learning algorithm for ELM includes three steps, which are as follows.Firstly, it determines the number of neurons of the hidden layer and sets the connection weight values  and the offset values  of the hidden layer neurons.Secondly, it selects an infinite differentiable function as an activation function of hidden layer neurons and computes the output array of the hidden layer .Finally, it calculates the weight value  of an output layer.Through the above steps, it randomly sets the weight and offset values of the hidden layer of ELM.It can also get the unique solution of output layer weights, and if the number of hidden layer nodes is enough, it can theoretically approach any continuous function.
The normally extreme learning machine produces the input weight values array and hidden layer offsets randomly, which will result in network instability and influence the forecast effect.Therefore, learning how to select the input weight values array in addition to the hidden layer offsets is key for improving ELM performance.Yu et al. proposed a method called Delta Test-ELM (DT-ELM) which can operate in an incremental way to create less complex ELM structures and determine the number of hidden nodes automatically [28].Shrivastava and Panigrahi proposed a hybrid wavelet-ELM based short term price forecasting for electricity markets [29].Focused on ELM performance problem, the paper introduced a global optimization intelligence algorithm, DE algorithm, to optimize the input weight values and hidden layer offsets in an effort to improve the algorithm.First, input weight values and hidden layer offsets are set as the optimizing variables of DE, and then the training forecast error is set as the fitness value of DE.Because the threedimensional numerical tests require too much time to perform, an orthogonal design method and a uniform design method are adopted for selecting computing schemes.They guarantee the accuracy of samples and also decrease the number of numerical calculation times.
The steps of the improved extreme learning machine are as follows.
Step 1.For each analysis task, training data sets are constructed to correspond to the geomechanical parameters and displacement at key points.Numerical analysis is used to calculate the data set for every set of orthogonal experimental schemes.To improve the generation performance of ELM, a testing sample set is selected from the data set and is used to assess the applicability of ELM.
Step 2. It sets the DE parameters, including population scale, evolution generation, cross factor, and magnification factor, and it randomly produces the first generation.The population scale is generally 20∼100, and each individual corresponds to the weight values and hidden layer offsets, training and obtaining the output weight values of ELM, and then it acquires the ELM structure.
Step 3.According to formulas (7), (8), and (9), each one forecasts and tests the trained ELM by test samples in order to set the maximal relative error as the fitness value of DE.The applicability of the model is measured in terms of fitness: where   and    are the estimated displacement of the tentative ELM and the calculated key point rock mass is at the th testing sample.The test number is  = 1, 2, . . ., .
Step 4. It carries out the variation, cross, and selection operations, and it iterates circularly until they meet the ending condition (maximal iteration).If the fitness is accepted, the ELM training procedure ends, outputting the input layer weight values and hidden layer offsets.If not, new population is produced according to (3), ( 4), (5), and (6), and the process returns to Step 3. (1), executing the DE iteration again.Then the geomechanical parameters are derived.
The study adopted ELM to express the relation between input and output data, instead of the numerical simulation, so as to save the computing time.The training data for extreme learning machine are required in an adequate and representative way.If the parameters are selected arbitrarily, too much numerical calculation and computing cost will also be induced.Accordingly, the following measures are considered: (1) orthogonal experimental and uniform designs are the methods for utilizing orthogonal table or uniform table to arrange and analyze the multiple factors test scientifically.Their main advantage is that they can select representative few test schemes among a lot of schemes.Adopting them to design parameters combination and constructing samples can improve test schemes and save test time, thus guaranteeing the samples representativeness.(2) Preliminary sensitivity means the degree of the factors influencing the system appraising indices.If the mechanics parameters are multiple, the preliminary sensitivity analysis should be applied to find out the most critical parameters affecting model characters which need back analysis so as to avoid unnecessary parameters, thereby increasing the calculation efficiency.

Engineering Introduction.
The Fusong tunnel is located at the Changbai Mountain of Jingyu County, Jilin Province of China.The exposed strata in the region are primarily kataclastics of fluviolacustrine of Jura, Andesite, and basalt of Tertiary and Quaternary systems.Moreover, the small-scale crack is not very concerned, and the degree of destruction is not serious.The underground water of the tunnel region is mainly phreatic water of loosening accumulative formation and bedrock fissure water.Additionally, the tunnel is designed as separate double caverns with a distance of 13∼ 35 m.The tunnel section has the height of 12 m and the width of 7.6 m.The left line has a length of 1625 m, starting at ZK275 + 170 and ending at ZK276 + 795.The longitudinal profile of the left line and the study region is shown in Figure 3.
According to the demand of highway tunnel construction norm of China, considering the surrounding rock and structure characteristics, the monitoring items were adopted, with an inclusion of top arc subsidence and convergence.And the automatic collection system including lining strain, rock pressure, and rock inner displacements was also monitored between the study regions.From November 27 to January 19, 2013, the field monitoring data curves and schematic diagrams were shown in Figures 4 and 5.

Three-Dimensional Numerical Simulation and Parameters
Sensibility Analysis.In order to carry out the parameters back analysis, firstly we constructed a three-dimensional numerical model for the left line tunnel of the engineering.The computing region is adopted with 207 meters in length, 70 meters in width, and 90 meters in height.The strata are cutin tuff, lime mudstone, and calcic silty mudstone.The coordinates of the tunnel midpoint are (0, 0, 0).-axis is the tunnel excavation direction, -axis is the vertical direction, and -axis is the horizontal radial direction.The model includes 370924 elements and 379464 nodes.
The strata adopted solid element and rock bolt; long pipe roof adopted cable structural element; shotcrete adopted shell element; and plastic constitutive model adapted to the Mohr-Coulomb criterion.According to actual monitoring, we set up a monitoring section every 20 meters.The threedimensional numerical model and typical section are shown in Figure 6.The selected parameters include elastic modulus, Poisson ratio, cohesion, and internal friction angle of lime mudstone ( 1 ,  1 ,  1 ,  1 ).In addition, the elastic modulus and Poisson ratio of cutin tuff ( 2 and  2 ) are study parameters.Five levels were adopted in the scope, and 25 orthogonal schemes of L25(6 5 ) were constructed.The results of orthogonal scheme numerical calculations are shown in Table 1.
In the back analysis process, most mechanics parameters are related to the surrounding rock stability and surrounding rock displacement, with different parameters having different effects on the surrounding rock deformation and destruction.The back analysis precision will decrease with the parameters increasing.So we carried out the parameters sensibility analysis and obtained the regulation of how the parameters affect the monitoring data.Based on the orthogonal scheme numerical calculations results (Table 1), through the extreme difference analysis, we know how the rock mechanics parameters affect stress and displacement of surrounding rock; the relations between parameters and difference ranges are, respectively, shown in Figures 7 and 8. From Figures 7 and 8, we obtained the parameters affecting the stress which in order are  1 >  1 >  1 >  2 >  2 >  1 and obtained the parameters with the effects on the displacement which in order are  1 >  1 >  1 >  2 >  1 >  2 .Combining the above results, elastic modulus  1 , cohesion  1 , and inner friction angle  1 are selected for succeeding analysis, with data fitting and back analysis based on DE-ELM.
Considering the forward numerical simulation costing too much time, using the DE-ELM algorithm to fit the relation between mechanics parameters and monitoring information, with ELM trained by the samples provided by Table 1, acquired the mapping expression of   = ().We tested the ELM regression model forecast effect by 5 uniform samples.The training time and mean relative error of forecast of the top arc subsidence compared with other similar methods named GA-SVM (genetic algorithm-support vector machine) and GA-BP (genetic algorithm-back propagation neural network) are listed in Table 2.
It is shown that DE-ELM algorithm has the best calculation speed and precision.In the ELM training process, the multiple inputting items and multiple outputting items are normalized according to the following formula: where   ,  min , and  max are, respectively, the th sample of the calculation value, the minimum value of samples, and the    factor  is 0.7, and the cross factor CR is 0.8.We executed the difference evolution algorithm, obtained the mechanics parameters  1 = 1.635GPa,  1 = 0.553 MPa, and  1 = 36.84∘ , and input the above back analyzed parameters into the numerical model.The calculated and field-monitored data are shown in Table 3.
It is shown that the displacement-stress coupling back analysis can get higher precision than the displacement back analysis; the feedback analysis relative errors of , , and  displacements changed from 8.608%, 10.313%, and 8.897% to 4.316%, 4.479, and 2.719% accordingly.

DE Parameters
Affecting the Optimizing Property.The inversion iterative curves of different strategies for the ELM optimizing are shown in Figure 9. DE/rand/1/bin and DE/rand/2/bin are selected; both of the two difference strategies will converge, but DE/rand/1/bin converges faster.The iterative charge is a more stable cable, and the convergence effect is better.Determining which difference strategy is suitable involves selecting the strategy that saves time, improves convergence speed, and also fits the method into the local optimal necessary conditions.Two important parameters are variation factor  and cross factor CR in DE.The DE/rand/1/bin difference strategy was chosen; in Figure 10, the iterative curve of the variation factor  = 0.9, CR = 0.5∼0.8,cross factor CR = 0.9, and  = 0.5∼0.8 is illustrated.The parameters selection affects evolution speed obviously, and the evolution stagnancy and premature phenomena will occur if inappropriate parameters are selected.By adjusting  and CR, or increasing population scale, the premature phenomena can be controlled.

Activation Function and Hidden Layer Nodes Affecting
ELM.We adopted the DE-LEM method for training 25 orthogonal samples and five forecasted test samples.The RMSE corresponding to different activation function is shown in Figure 11.It is shown that the forecast error corresponding to the sigmoid activation function is the minimum, and the forecast error corresponding to the hardlim activation function is the maximum.

Conclusion
The paper conducted a research on the surrounding rock back analysis method based on the DE-ELM algorithm for highway tunnel construction and arrived at the following conclusions: (1) We introduced a new single hidden layer feed forward neural networks machine learning method, ELM, by combining it with the DE algorithm.The input weight values and hidden layer offsets were optimized, and the forecast precision was improved.We then constructed the surrounding rock identification method.The method took advantage of DE's global optimization property and the simple structure of ELM.The functions of pattern recognition and data fitting of ELM could be used for rock classification and parameter identification, respectively.
(2) The case study stated that the DE-ELM method achieved satisfactory identification with precision and high calculation speed.The activation function affected the forecast precision.The forecast error corresponding to the sigmoid activation function is the minimum, and the forecast error corresponding to hardlim activation function is the maximum; the difference between them is less than 0.1.The DE-ELM calculation results that corresponded to different hidden layer nodes were compared; 150-200 hidden layer nodes were recommended.
(3) The DE-ELM has been used in the construction of the Fusong tunnel where it ascertained the surrounding rock parameters for the region near YK275 + 370;

Figure 1 :
Figure 1: The forward and back analysis of geotechnical engineering.

Figure 2 :
Figure 2: The generating process of V ,+1 in two-dimension solution space.

Step 5 .Horizontal scale 1 : 2000 Figure 3 :
Figure 3: The longitudinal profile of the left line and the study area.

Figure 4 :
Figure 4: The field monitoring data curves of displacement measuring point and schematic diagram.

Figure 5 :
Figure 5: The field monitoring data curves of earth pressure cell and schematic diagram.

Figure 6 :
Figure 6: The 3D numerical model and section.
is the output value of the th sample,   = [ 1 ,  2 , . . .,   ] T expresses the connection weight values from the input layer to the hidden layer,   = [ 1 ,  2 , . . .,   ] T expresses the offset value,   = [ 1 ,  2 , . . .,   ] T expresses the connection weight values from hidden layer to output layer, and () is the active function that can adopt sigmoid function, Sine function, RBF function, and so forth.If the networks approach training samples with zero error, then

Table 1 :
The results of orthogonal scheme numerical calculation.

Table 2 :
Comparison between test and forecast samples of top arc subsidence.

Table 3 :
The comparison between unitary displacement feedback analysis and joint feedback analysis.