Abstract
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 algorithm proposed. An extreme learning machine (ELM) is employed with 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 by DE optimization arithmetic taking into consideration the monitoring 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.
1. 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–8]. Gioda et al. studied and compared the back analysis procedures for the interpretation of field measurements in geomechanics [9–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 three-dimensional 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 multi-information (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.
2. 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 () 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, is monitoring displacements (or pressure, stress) of tunnel-surrounding rock, is calculation displacements of tunnel-surrounding rock, and is the number of monitoring points.
2.2. 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 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 , . 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. Considerwhere rand is a random number within and 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 , , and are a random integer from , which is different from the others. 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 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 cross based on the following rules are to generate a trial vector : wherewhere is the component of the vector corresponding to random number, is the crossover factor, and is an integer randomly selected in .
(4) Selection. The greedy search method used to select the new solution vector 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.
2.3. 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 (), the th sample output value of SLFN with hidden layer nodes and a hidden layer activation function of is as follows:where is the output value of the th sample, expresses the connection weight values from the input layer to the hidden layer, expresses the offset value, 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 It can be simply expressed as wherewhere , , 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 . 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 three-dimensional 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 .
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.
Step 5. Set formula (1) as the fitness function of DE; adopt the successfully trained ELM of Step 4 as of formula (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.
3. Case Study
3.1. 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.
3.2. 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 (). -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 three-dimensional 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 (). In addition, the elastic modulus and Poisson ratio of cutin tuff ( and ) are study parameters. Five levels were adopted in the scope, and 25 orthogonal schemes of L25(65) were constructed. The results of orthogonal scheme numerical calculations are shown in Table 1.
(a) Section of geological body
(b) The 3D numerical model
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 and obtained the parameters with the effects on the displacement which in order are . Combining the above results, elastic modulus , cohesion , and inner friction angle 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 , , and are, respectively, the th sample of the calculation value, the minimum value of samples, and the maximum value of samples and is the normalization of the th samples, .
Input the derived ELM model and the field monitoring data mm, mm, mm, KPa, KPa, and KPa; set the DE parameters optimizing variables as 3; the population scale is 50, the iteration number is 600, the variation factor is 0.7, and the cross factor CR is 0.8. We executed the difference evolution algorithm, obtained the mechanics parameters GPa, MPa, and , 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.
4. Discussion
4.1. 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 , 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.
4.2. 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.
Figure 11 also shows the forecast RMSE variation according to the increasing number of hidden layer nodes. The RMSE corresponding to 30 hidden layer nodes’ number is larger than 0.5, whereas when there are 150 hidden layer nodes, the corresponding RMSE approaches 0.2. However, as the hidden layer nodes increase from 150 to 270, the RMSE does not change obviously. From the results, in the similar back analysis calculation by DE-LEM, 150–200 hidden layer nodes are recommended.
5. 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; what is stated is that the rock is IV class, and the identified surrounding rock parameters are = 3.496 GPa, = 0.315, and = 0.300 MPa. The adjacent region construction scheme changed from the Two Side-Wall Pilot Tunnel Method to Benching Tunneling Construction Method. The use of the adjusted construction scheme realized the expected target, demonstrating that the proposed back analysis method is reasonable.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors deeply appreciate support from the Fundamental Research Funds for the Central Universities (2011JC012, 2012TD015) and the National Natural Science Foundation (51079010).