Tunnel surrounding rock parameter intelligent inversion analysis method based on XGBoost optimization algorithm

Abstract

This invention belongs to the field of tunnel engineering stability analysis, specifically disclosing an intelligent inversion analysis method for tunnel surrounding rock parameters based on the XGBoost optimization algorithm, including the following steps: S1: Establish a numerical simulation calculation model to obtain a sample library for surrounding rock parameter inversion; S2: Perform correlation and sensitivity analysis on the displacement and parameter data in the sample library, evaluate the feasibility of each parameter as the parameter to be inverted, and quantitatively evaluate the rationality of displacement feature combinations; S3: Use the CART algorithm to determine the parameters to be inverted and perform displacement feature combination screening; S4: Use the CART algorithm as the base learner to establish an XGBoost ensemble algorithm model for intelligent inversion analysis of tunnel surrounding rock parameters; S5: Apply the single variable control method and Bayesian optimization method to optimize the XGBoost algorithm; S6: Input the displacement features into the XGBoost ensemble algorithm model to obtain the predicted surrounding rock parameter values. This invention achieves optimization of the XGBoost ensemble algorithm model through hyperparameter optimization, resulting in high model stability and prediction accuracy.

Description

Technical Field

[0001] This invention belongs to the field of tunnel stability analysis, specifically a method for intelligent inversion analysis of tunnel surrounding rock parameters based on the XGBoost optimization algorithm. Background Technology

[0002] Numerical simulation is a commonly used method for tunnel stability analysis, as it can easily obtain the required results such as tunnel surrounding rock stress, displacement, and stability. However, the accuracy of the calculation results from numerical simulation models is largely dependent on the input values ​​of the surrounding rock parameters. In actual tunnel engineering surveys and construction, due to factors such as project cost and construction period, it is difficult to obtain tunnel surrounding rock parameters accurately and in a timely manner through laboratory tests and in-situ tests.

[0003] The tunnel surrounding rock inversion analysis method, which uses displacement values ​​to invert surrounding rock parameters, was proposed in the 1970s. After decades of continuous development and practice, numerous research findings have proven the feasibility and effectiveness of displacement-based surrounding rock parameter inversion analysis.

[0004] With the development of mathematics and computer science, the method of using machine learning algorithms to process complex numerical variable data and explore the mapping relationships between variables has begun to be applied in various fields of production and life. From the late 1990s to the present, using numerical simulation analysis and various machine learning algorithms for tunnel surrounding rock parameter inversion analysis has become the mainstream method in the field of tunnel stability analysis.

[0005] Typically, the technical process for parameter inversion using numerical simulation and machine learning algorithms is as follows:

[0006] 1. Establish a three-dimensional numerical simulation model for tunnel excavation and support of the tunnel section to be analyzed;

[0007] 2. Determine the design combination of mechanical parameters as needed (e.g., use orthogonal test schemes to design combinations of parameters such as elastic modulus, Poisson's ratio, unit weight, cohesion, and internal friction angle), substitute the parameter combination into the numerical simulation calculation model to calculate and extract multiple sets of displacement characteristic values, and establish a sample library of surrounding rock parameters-displacement;

[0008] 3. Determine the parameters to be inverted and perform displacement feature screening;

[0009] 4. Based on machine learning algorithms, establish the mapping relationship between surrounding rock displacement characteristics and parameters to be inverted. Divide the sample database into training and test sets. Use the training set data to train the machine learning algorithm model, and use the test set data to verify the accuracy and usability of the model.

[0010] 5. By substituting the measured surrounding rock displacement into the trained model, the predicted surrounding rock parameter values ​​can be obtained.

[0011] As can be seen from the above process, the key issue in this technical work lies in the machine learning algorithm employed. Ensuring the reliability of the sample database while improving the prediction accuracy and stability of the parameter inversion model is the focus of this technical work.

[0012] Previous research has extensively explored the use of displacement values ​​to invert tunnel surrounding rock parameters, employing various machine learning algorithms such as artificial neural networks, immune Gaussian processes, differential evolution, and support vector machines. However, previous studies have generally relied on single algorithms, which suffer from limited prediction accuracy and cannot fully leverage the algorithm's potential. While existing research has also used the XGBoost ensemble algorithm for similar analyses, it lacks hyperparameter optimization, resulting in low model stability and prediction accuracy. Furthermore, current research has not evaluated the feasibility of using various surrounding rock parameters as inversion parameters or the rationality of combining displacement features. Summary of the Invention

[0013] The purpose of this invention is to address the low prediction accuracy issues in current research on surrounding rock parameter inversion analysis using tunnel displacement values, and to propose using the XGBoost ensemble algorithm to improve the model's prediction accuracy. In parameter inversion based on machine learning algorithms, not all surrounding rock parameters are suitable as inversion parameters; therefore, this invention evaluates the feasibility of using each parameter as an inversion parameter. To address the problem of poor stability of the inversion model due to high correlation of displacement feature values, this invention performs displacement feature screening. Regarding the optimization of the ensemble algorithm itself, this invention applies the controlled single variable method and Bayesian optimization method to optimize the hyperparameters of the XGBoost ensemble algorithm, thereby improving the model's stability and predictive ability.

[0014] This invention discloses an intelligent inversion analysis method for tunnel surrounding rock parameters based on the XGBoost optimization algorithm, comprising the following steps:

[0015] S1: Establish a three-dimensional finite element numerical simulation model for tunnel excavation and support, and obtain a sample library for inverting surrounding rock parameters;

[0016] S2: Perform feature correlation and sensitivity analysis on displacement and parameter data in the sample library, evaluate the feasibility of each parameter as the parameter to be inverted, and evaluate the rationality of displacement feature combination.

[0017] S3: Use the CART algorithm to perform parameter inversion analysis and select displacement features;

[0018] S4: Using the CART algorithm as the base learner, establish an XGBoost integrated algorithm model for intelligent inversion analysis of tunnel surrounding rock parameters;

[0019] S5: Apply the single-variable control method and Bayesian optimization method to perform hyperparameter optimization, and obtain the XGBoost integrated algorithm model for intelligent inversion analysis of optimized tunnel surrounding rock parameters;

[0020] S6: Input the displacement features into the trained XGBoost model for intelligent inversion analysis of surrounding rock parameters to obtain the predicted values ​​of surrounding rock parameters.

[0021] In a preferred embodiment of the present invention, in S1, an orthogonal experimental design is performed on the surrounding rock parameters of the rock stratum where the tunnel is located to obtain multiple sets of surrounding rock parameter combinations; the multiple sets of surrounding rock parameters are input into the numerical simulation calculation model, and the displacement values ​​of the surrounding rock deformation stabilization stage are extracted to obtain a sample library for surrounding rock parameter inversion.

[0022] In a preferred embodiment of the present invention, in S2, the feature correlation analysis is calculated using the Pearson correlation coefficient. The Pearson correlation coefficient between two variables X and Y is the quotient of the covariance and standard deviation between the two variables. The closer the absolute value of the Pearson correlation coefficient is to 1, the stronger the linear relationship between the two features. The formula for the Pearson correlation coefficient is:

[0023]

[0024] In a preferred embodiment of the present invention, in S2, the average value of the parameters of the rock strata where the tunnel is located, 90% of the average value of the parameters, and 110% of the average value of the parameters are taken and substituted into the numerical simulation calculation model to calculate the corresponding displacement. By comparing the displacement changes when different parameters change by the same proportion, the parameter sensitivity of the surrounding rock is determined.

[0025] As a preferred embodiment of the present invention, in S3, the CART algorithm is a binary decision tree model, and its specific calculation process is as follows:

[0026] (1) Select the segmentation variable j, where the segmentation variable is the characteristic value of each displacement;

[0027] (2) Select the dividing point s according to the dividing principle;

[0028] (3) For any indivisible variable and any indivisible point s i The training sample x is split into two sample subsets R. 1i and R 2i ;

[0029] Among them, R 1i ={x j |x j ≤s i}, j=1~n;R 2i ={x j |x j >s i}, j = 1 to n; i is the number of the split point;

[0030] (4) Let c1 and c2 be R respectively. 1i and R 2i The average value of the target value y in the middle sample:

[0031]

[0032]

[0033] (5) Based on the principle of minimizing the squared error, the optimal splitting variables and splitting points are selected:

[0034]

[0035] (6) Based on the selected splitting variable j and splitting point s, the training samples are divided into two sample subsets, R1 and R2.

[0036] (7) Repeat steps (1) to (6) in the sample subsets R1 and R2 to further divide the subsets R1 and R2 into smaller subsets;

[0037] (8) Repeat steps (1) to (7) to further divide the subset until the termination condition is met (such as reaching the maximum depth of the tree, the number of samples corresponding to the leaf node reaching the minimum number of samples, etc.).

[0038] (9) Finally, through the above division, the training sample input space is divided into R1, R2, R3, ... R m With m leaf nodes, the production decision tree is formed.

[0039]

[0040]

[0041]

[0042] In a preferred embodiment of the present invention, in S4, using the CART algorithm as the base learner, the XGBoost ensemble algorithm model is constructed as follows:

[0043] (1) The CART algorithm is iterated. Each iteration fits the residual of the CART decision tree obtained in the previous iteration to minimize the residual. The residual calculation formula is:

[0044] r = yf t-1 (x)

[0045] Where r is the residual, y is the true value of the sample output, t is the number of iterations, x is the input variable of the sample, and f t-1 (x) represents the predicted value of the CART decision tree model at the (t-1)th iteration;

[0046] (2) Perform the 0th iteration:

[0047]

[0048] Where f0(x) is the predicted value in the 0th iteration, L is the objective function, and y i Here, c represents the true value of the sample output variable, and c represents the parameter of the objective function. In this patent, c is taken as the average value of the sample output variable.

[0049] (3) Perform the first iteration: Use the difference between the true value of the sample output variable and the average value c of the sample output variable as the sample output value y. 1i Based on the specified objective function, the CART algorithm is applied to construct the first CART decision tree model and obtain the predicted value f1(x) of the sample output variable.

[0050] (4) Perform the second iteration: using the original output value y from the first iteration sample. 1i The difference between the predicted value f1(x) of the output variable of the first iteration sample and the predicted value f1(x) is taken as the output value y of the sample. 2i Based on the specified objective function, the CART algorithm is applied to construct the second CART decision tree model and obtain the predicted value f2(x) of the sample output variable.

[0051] (5) Following the same method as the second iteration, construct the 3rd, 4th, ..., Mth CART decision tree models to obtain the predicted values ​​f3(x), f4(x), ..., f5(x) of the sample output variables. M (x), until the residual meets the requirements;

[0052] (6) Based on the above M CART decision tree models, calculate the output value of the predicted sample according to the following formula:

[0053]

[0054] In a preferred embodiment of the present invention, in S4, the objective function of the XGBoost algorithm is obtained by adding a regularization term to the conventional squared error loss function:

[0055]

[0056] Where T is the number of leaf nodes; γ and λ are the regularization coefficients for leaf nodes and leaf node weights, respectively. It is the sum of the first derivatives of all samples under leaf node j; This is the first derivative of the loss function with respect to the predicted values ​​of the sample output variables; It is the sum of the second derivatives of the j samples under the leaf node; I is the second derivative of the loss function with respect to the predicted values ​​of the sample output variables; j ={i|q(x i Let ) = j} be the sample set under leaf node j; It is the prediction value for the (t-1)th round given after the accumulation of the learner in the first t-1 rounds; y i is the original value of the sample output variable; i is the sample number.

[0057] In a preferred embodiment of the present invention, in S5, the range of hyperparameter values ​​is narrowed by controlling a single variable; within the range narrowed by controlling a single variable, the optimal value of the hyperparameter is further determined by Bayesian optimization.

[0058] The specific process of the single-variable method is as follows:

[0059] (1) For each hyperparameter, the hyperparameter values ​​are discretized to obtain a certain number of equal points;

[0060] (2) Substitute the discretized hyperparameter values ​​into the XGBoost ensemble algorithm model to obtain the coefficient of determination and generalization error of each model;

[0061] (3) Based on the variation law of the coefficient of determination and the generalization error, obtain the optimal range of values ​​for each hyperparameter.

[0062] Within the narrowed hyperparameter range, the process of finding hyperparameters using Bayesian optimization is as follows: (1) Establish a hyperparameter sample as x = [x1, x2, x3, ..., x n The objective function of the XGBoost ensemble algorithm model takes the value y = [y1, y2, y3, ..., y]. n ], y follows a multidimensional normal distribution;

[0063] (2) Establish the kernel function for the Gaussian process:

[0064]

[0065] in, and l are the parameters of the kernel function;

[0066] (3) Based on the kernel function, establish the kernel vector:

[0067]

[0068] (4) Establish the likelihood function of y:

[0069]

[0070] Where μ is the mean of x;

[0071] The optimal kernel matrix is ​​obtained by minimizing the likelihood function using gradient descent and calculating the parameters of the kernel function.

[0072] (5) Based on the optimal kernel matrix, establish a Gaussian process between the posterior probability and prior probability of y and the input variables:

[0073]

[0074] Where y * K is the value of the output variable obtained in the next iteration. * This is the kernel vector for the next iteration;

[0075] The Gaussian process described above establishes the functional relationship between the objective function value and the hyperparameters of the model.

[0076] (6) Using the above Gaussian process as a probabilistic surrogate model, establish POI as the learned function:

[0077]

[0078] Where f(X) is the objective function value of X, i.e., the Gaussian process derived above; f(X) + ) represents the current optimal objective function value, i.e., the posterior distribution of f(X); μ(x) and σ(x) are the mean and variance of the objective function obtained by the Gaussian process; ξ is the trade-off coefficient, which controls whether the optimization direction of the hyperparameters is biased towards "development" or "search"; based on the above POI learned function, the Monte Carlo method is applied to search for the above optimal hyperparameters.

[0079] (7) Substitute the optimal hyperparameters searched based on the learned function into the probabilistic proxy model to establish a new Gaussian process;

[0080] (8) Using the new Gaussian process as the probabilistic surrogate model, establish a new POI and apply the Monte Carlo method to perform a new round of hyperparameter optimization.

[0081] (9) Repeat the above operation until the iteration stopping condition is met to obtain the final optimal hyperparameter combination.

[0082] Explanation: POI is an abbreviation for Probability of Improvement.

[0083] In summary, the outstanding technical advantages of this invention are:

[0084] 1) Conducting feature correlation analysis and parameter sensitivity analysis provides support for determining the parameters to be inverted, because not all surrounding rock parameters can be inverted using displacement.

[0085] 2) Filtering displacement features can improve the stability of the model, minimize the adverse effects of high correlation between displacement features on the model, and also consider the prediction accuracy to select the optimal combination of displacement features.

[0086] 3) As a nonlinear algorithm, the CART algorithm can effectively avoid the adverse effects of the high correlation of displacement characteristics on the model in linear model algorithms. Compared with linear algorithms, the CART algorithm has higher prediction accuracy in parameter inversion applications and has better applicability.

[0087] 4) The XGBoost ensemble algorithm model for intelligent analysis of tunnel surrounding rock parameters inversion based on the CART algorithm can give full play to the iterative optimization role of the ensemble algorithm. The XGBoost ensemble algorithm focuses on the nodes that are misclassified in each iteration of the CART algorithm and adjusts the weight values ​​of the nodes, which can effectively improve the prediction accuracy of the base learner CART algorithm.

[0088] 5) Combining the single-variable control method and Bayesian optimization to optimize the hyperparameters of the XGBoost ensemble algorithm model can significantly improve the model's prediction accuracy and stability, as well as its fitting effect. Compared with traditional single-method approaches, such as grid search, random search, and Bayesian optimization, the combined optimization method is more efficient and accurate. Attached Figure Description

[0089] Figure 1 The flowchart is a process for the intelligent inversion analysis method of tunnel surrounding rock parameters based on the XGBoost optimization algorithm of this invention.

[0090] Figure 2 This is a schematic diagram of the rock strata division in a three-dimensional model of tunnel excavation in an embodiment of the present invention;

[0091] Figure 3 This is a schematic diagram of the detailed grid dimensions for tunnel excavation in an embodiment of the present invention;

[0092] Figure 4 This is a stress diagram of the surrounding rock of the tunnel after excavation in an embodiment of the present invention;

[0093] Figure 5 This is a strain cloud diagram of the tunnel surrounding rock after excavation in an embodiment of the present invention;

[0094] Figure 6 This is a schematic diagram of the displacement extraction points in an embodiment of the present invention;

[0095] Figure 7This is a flowchart of the CART algorithm in the intelligent inversion analysis method for tunnel surrounding rock parameters based on the XGBoost optimization algorithm of this invention;

[0096] Figure 8 This is a flowchart of the XGBoost integrated algorithm in the intelligent inversion analysis method for tunnel surrounding rock parameters based on the XGBoost optimization algorithm of this invention;

[0097] Figure 9 This is a graph showing the coefficient of determination scores of the hyperparameter n_estimators when it takes different values ​​in an embodiment of the present invention.

[0098] Figure 10 This is a curve showing the generalization error score when the hyperparameter n_estimators takes different values ​​in an embodiment of the present invention.

[0099] Figure 11 This is a bar chart comparing the predicted and actual values ​​of Poisson's ratio in an embodiment of the present invention. Detailed Implementation

[0100] The present invention will now be described in detail with reference to the accompanying drawings.

[0101] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0102] This embodiment presents an intelligent inversion analysis method for tunnel surrounding rock parameters based on the XGBoost optimization algorithm, as shown in the attached figure. Figure 1 As shown, the process includes: S1: Using finite element software, a three-dimensional finite element numerical simulation model for tunnel excavation is established to obtain a sample library for the inversion of surrounding rock parameters; S2: Feature correlation and sensitivity analysis is performed on the displacement and parameter data in the sample library to evaluate the feasibility of each parameter as a parameter to be inverted and the rationality of the displacement feature combination; S3: Parameter inversion analysis is performed using the CART algorithm to screen displacement features; S4: Using the CART algorithm as the base learner, an XGBoost ensemble algorithm model for intelligent inversion analysis of tunnel surrounding rock parameters is established; S5: Hyperparameter optimization is performed using the controlled single variable method and Bayesian optimization method to obtain an optimized XGBoost ensemble algorithm model for intelligent inversion analysis of tunnel surrounding rock parameters; S6: The displacement features are input into the trained XGBoost model for intelligent inversion analysis of surrounding rock parameters to obtain the predicted surrounding rock parameter values.

[0103] The following specific examples will provide a more detailed description of the present invention:

[0104] The first step involves modeling the right-line tunnel (a separate mainline tunnel from K4+173 to K4+205) of a certain tunnel project, studying the surrounding rock displacement and parameters after excavation and support stabilization at section K4+173. This section of the tunnel has a depth of 154 meters, with high-strength surrounding rock classified as Class III. Drill-and-blast method is being used for construction, and the right-line tunnel has been advanced to a sufficiently far position. The distance between the left and right tunnels is considerable, meeting the spacing requirements for a separated mainline tunnel, and their mutual influence can be ignored. A full-face excavation method without support is employed to numerically simulate the tunnel excavation and support construction stages.

[0105] The displacement data in the sample library is extracted based on the size of the specific excavation section and the support conditions. It generally includes multiple sets of displacement characteristics such as crown settlement, arch bottom heave, and horizontal convergence. The parameters in the sample library are the input parameters when using the Mohr-Coulomb constitutive model for the soil and rock mass, including five parameters: elastic modulus, Poisson's ratio, unit weight, cohesion, and internal friction angle. The time required for the surrounding rock to enter the stable stage is determined based on the measured displacement time-series curve; that is, the extracted stable displacement value is uniformly determined as the stable displacement value on the 32nd day after tunnel excavation. The tunnel is constructed in one step per day, with each step excavating 1 meter. Therefore, the extension distance of the tunnel excavation 3D model in the depth direction is 32 meters. In the width and height directions, the extension range must be approximately five times the excavation size. In this example, the horizontal length of the constructed tunnel excavation 3D model is 199.5 meters, and the vertical height is 218.8 meters. The 3D model uses a hexahedral mesh, with the mesh size for the tunnel excavation section and its surrounding rock strata controlled at 0.6 meters, the mesh size for adjacent rock strata controlled at 0.8 meters, and the mesh size for the remaining rock strata controlled at 1.2 meters. The numerical simulation model for the tunnel excavation was obtained using finite element method (FEM) software, including but not limited to Midas GTS NX, FLAC 3D, and ABAQUS. The constitutive relationship of the soil and rock mass in the numerical simulation model adopts the Mohr-Coulomb constitutive model. A schematic diagram of the rock strata division in the constructed 3D model is attached. Figure 2 As shown in the attached diagram, the detailed grid dimensions for tunnel excavation are as follows. Figure 3 As shown.

[0106] The preliminary geological survey report indicates that the tunnel section has six rock strata, from top to bottom: strongly weathered sandstone, moderately weathered sandstone, black shale, gray sandstone, black shale, and gray sandstone. During modeling, the surrounding rock parameters for the first five strata were taken from the average values ​​of those parameters in the geological survey report. The sixth stratum, gray sandstone, is the stratum where the tunnel excavation will take place, and its parameters were input according to orthogonal design. The maximum and minimum values ​​of these parameters were determined based on the survey report and actual conditions to define the range of values ​​for each parameter. The parameter values ​​for each stratum are shown in Table 1.

[0107] Table 1

[0108]

[0109] Step 2: Orthogonal experimental design was conducted on five parameters of the gray sandstone stratum where the tunnel excavation is located. Based on the value ranges of each parameter in the stratum determined in Step 1, an orthogonal design was performed. In this example, SPSS software was used, with six value levels for each parameter, for a total of five parameters. Therefore, the orthogonal design generated 49 combinations of surrounding rock parameters, as shown in Table 2, which presents some of the five-factor, six-level orthogonal experimental parameter combinations.

[0110] Table 2

[0111]

[0112] Determine the location of the displacement feature extraction points. Select the mean or median of the parameters of the surrounding rock strata where the tunnel is excavated, substitute them into the numerical simulation model, and observe the stress and strain contour maps after the tunnel excavation. The stress diagram is shown below. Figure 4 As shown, the strain contour plot is as follows Figure 5 As shown in the figure, displacement monitoring points are reasonably arranged in areas with larger stress-strain results. The displacements of the monitoring points are as follows: Figure 6 As shown.

[0113] The 49 sets of surrounding rock parameters designed by orthogonal experiment were successively substituted into the numerical simulation calculation model for calculation, and the displacement values ​​of the corresponding extraction points were extracted to form a parameter-displacement sample library for parameter inversion. Some parameter-displacement data are statistically summarized in Table 3.

[0114] Table 3

[0115]

[0116] Step 3: Perform feature correlation analysis on the parameter-displacement data in the database. Feature correlation analysis uses the Pearson correlation coefficient for calculation. The Pearson correlation coefficient between two variables X and Y is the quotient of their covariance and standard deviation. The closer the absolute value of the Pearson correlation coefficient is to 1, the stronger the linear relationship between the two features. The formula for the Pearson correlation coefficient is:

[0117]

[0118] The Pearson correlation coefficient scores among the variables are shown in Table 4 below.

[0119] Table 4

[0120]

[0121]

[0122] Step 4: Perform parameter sensitivity analysis. First, calculate the average value of each parameter provided in the geological survey report for the rock strata where the tunnel excavation is located (layer 6: gray sandstone). Then, calculate 0.9 times and 1.1 times for each parameter. Using the controlled single variable method, combine 0.9 times, 1 times, and 1.1 times for one parameter with 1 times for another parameter each time, obtaining three parameter combination forms each time. Substitute these parameter combinations into the numerical simulation model of tunnel excavation for calculation, extract displacement characteristic values, and then subtract the displacement value generated by 0.9 times the parameter from the displacement characteristic value generated by 1.1 times the parameter. Take the absolute value of the difference, divide the absolute value of the difference by the displacement characteristic value generated by 1 times the parameter, and obtain the percentage change of each displacement characteristic when each parameter changes by the same proportion. This percentage quantifies the magnitude of displacement generated when the parameter changes. Taking unit weight and elastic modulus as examples, the calculation process of parameter sensitivity is shown in Table 5, and the calculation results of parameter sensitivity are shown in Table 6.

[0123] Table 5

[0124]

[0125] Table 6

[0126]

[0127]

[0128] Step 5: Analysis based on the results of variable correlation and parameter sensitivity analysis: The results of variable correlation analysis show that the correlation scores between elastic modulus, Poisson's ratio, and displacement characteristics are high, indicating a strong correlation between these two parameters and displacement characteristics; the correlation scores between unit weight, cohesion, internal friction angle, and displacement characteristics are low, indicating a weak correlation between these three parameters and displacement. The results of parameter sensitivity analysis show that changes in elastic modulus and Poisson's ratio result in significant displacement changes, changes in unit weight result in very small displacement changes, and changes in cohesion and internal friction angle do not result in any displacement changes. Based on the combined results of parameter sensitivity and variable correlation analysis, the parameters to be inverted are preliminarily determined to be elastic modulus and Poisson's ratio.

[0129] The analysis of the correlation calculations revealed a high linear correlation between the independent variables, specifically the displacement characteristics. The correlation coefficient between convergence two and convergence three was 0.96, and the correlation coefficient between the crown drop and the base rise was 0.98. Given this high correlation, using a linear algorithm would result in poor model reliability and low confidence. Therefore, a nonlinear algorithm was chosen for the parameter inversion analysis.

[0130] Step 6: Use the CART algorithm to determine the parameters to be inverted and perform displacement feature combination screening. The CART algorithm flowchart is attached. Figure 7 As shown, the CART algorithm, as a nonlinear algorithm, can effectively avoid the adverse effects of linear correlation and multicollinearity between displacement features in linear algorithm models. CART regression is a recursive process of generating a binary tree, and the splitting process uses the criterion of minimizing the squared error. Its specific calculation process is as follows:

[0131] (1) Select the segmentation variable j, where the segmentation variable is the characteristic value of each displacement;

[0132] (2) Select the dividing point s according to the dividing principle;

[0133] (3) For any indivisible variable and any indivisible point s i The training sample x is split into two sample subsets R. 1i and R 2i ;

[0134] Among them, R 1i ={x j |x j ≤s i}, j=1~n;R 2i ={x j |x j >s i}, j = 1 to n; i is the number of the split point;

[0135] (4) Let c1 and c2 be R respectively. 1i and R 2i The average value of the target value y in the middle sample:

[0136]

[0137]

[0138] (5) Based on the principle of minimizing the squared error, the optimal splitting variables and splitting points are selected:

[0139]

[0140] (6) Based on the selected splitting variable j and splitting point s, the training samples are divided into two sample subsets, R1 and R2.

[0141] (7) Repeat steps (1) to (6) in the sample subsets R1 and R2 to further divide the subsets R1 and R2 into smaller subsets;

[0142] (8) Repeat steps (1) to (7) to further divide the subset until the termination condition is met (such as reaching the maximum depth of the tree, the number of samples corresponding to the leaf node reaching the minimum number of samples, etc.).

[0143] (9) Finally, through the above division, the training sample input space is divided into R1, R2, R3, ... R m With m leaf nodes, the production decision tree is formed.

[0144]

[0145]

[0146]

[0147] Table 7 shows the statistical results of parameter inversion using the CART algorithm. The coefficient of determination (r) is used to characterize the accuracy of the parameter inversion. 2 The calculation formula is as follows:

[0148]

[0149] Table 7

[0150]

[0151] As shown in Table 7, the CART algorithm has a low accuracy rate in inverting the elastic modulus, so the feasibility of using the CART algorithm to invert the elastic modulus is poor. However, the CART algorithm has a high accuracy rate in inverting Poisson's ratio, so the feasibility of establishing a prediction model for inverting Poisson's ratio is relatively high.

[0152] To improve the prediction accuracy of the CART algorithm model, the XGBoost ensemble algorithm was used to integrate the individual algorithms. However, the XGBoost ensemble algorithm requires the prediction accuracy of each individual algorithm to be higher than 0.5. Therefore, the parameter to be inverted was further determined to be Poisson's ratio.

[0153] To further determine the optimal combination of displacement features, the displacement features were screened. Specifically, five displacement feature combinations were enumerated: crown settlement, arch bottom heave, horizontal convergence one, horizontal convergence two, and horizontal convergence three, resulting in a total of 19 displacement parameter combinations. The calculation results of Poisson's ratio from these parameter combinations are shown in Table 8.

[0154] Table 8

[0155]

[0156] The calculation results of the Poisson's ratio inversion from each displacement feature combination show that the prediction accuracy is the highest when the feature combination is selected as all displacement features. Therefore, this displacement feature combination is the optimal displacement feature combination.

[0157] Step 7: Using the CART algorithm as the base learner, establish a tunnel surrounding rock parameter inversion model based on the XGBoost ensemble algorithm. The XGBoost ensemble algorithm flowchart is attached. Figure 8As shown, XGBoost is an ensemble algorithm that integrates base learner algorithms based on the Boosting concept. In each iteration, XGBoost adjusts the weights of the leaf nodes of the CART base learner, focusing on nodes that were misclassified in the previous iteration and modifying their weights. The final prediction result of XGBoost is the sum of the weights of all nodes.

[0158] The specific implementation process for constructing the XGBoost ensemble algorithm model is as follows:

[0159] (1) The CART algorithm is iterated. Each iteration fits the residual of the CART decision tree obtained in the previous iteration to minimize the residual. The residual calculation formula is:

[0160] r = yf t-1 (x)

[0161] Where r is the residual, y is the true value of the sample output, t is the number of iterations, x is the input variable of the sample, and f t-1 (x) represents the predicted value of the CART decision tree model at the (t-1)th iteration;

[0162] (2) Perform the 0th iteration:

[0163]

[0164] Where f0(x) is the predicted value in the 0th iteration, L is the objective function, and y i Here, c represents the true value of the sample output variable, and c represents the parameter of the objective function. In this patent, c is taken as the average value of the sample output variable.

[0165] (3) Perform the first iteration: Use the difference between the true value of the sample output variable and the average value c of the sample output variable as the sample output value y. 1i Based on the specified objective function, the CART algorithm is applied to construct the first CART decision tree model and obtain the predicted value f1(x) of the sample output variable.

[0166] (4) Perform the second iteration: using the original output value y from the first iteration sample. 1i The difference between the predicted value f1(x) of the output variable of the first iteration sample and the predicted value f1(x) is taken as the output value y of the sample. 2i Based on the specified objective function, the CART algorithm is applied to construct the second CART decision tree model and obtain the predicted value f2(x) of the sample output variable.

[0167] (5) Following the same method as the second iteration, construct the 3rd, 4th, ..., Mth CART decision tree models to obtain the predicted values ​​f3(x), f4(x), ..., f5(x) of the sample output variables. M(x), until the residual meets the requirements;

[0168] (6) Based on the above M CART decision tree models, calculate the output value of the predicted sample according to the following formula:

[0169]

[0170] The objective function of the XGBoost algorithm is obtained by adding a regularization term to the conventional squared error loss function:

[0171]

[0172] Where T is the number of leaf nodes; γ and λ are the regularization coefficients for leaf nodes and leaf node weights, respectively. It is the sum of the first derivatives of all samples under leaf node j; This is the first derivative of the loss function with respect to the predicted values ​​of the sample output variables; It is the sum of the second derivatives of the j samples under the leaf node; I is the second derivative of the loss function with respect to the predicted values ​​of the sample output variables; j ={i|q(x i Let ) = j} be the sample set under leaf node j; It is the prediction value for the (t-1)th round given after the accumulation of the learner in the first t-1 rounds; y i is the original value of the sample output variable; i is the sample number.

[0173] Table 9 is a statistical table evaluating the performance of the XGBoost ensemble algorithm in retrieving Poisson's ratio.

[0174] Table 9

[0175]

[0176] The XGBoost ensemble algorithm has numerous hyperparameters. When faced with a specific problem, the characteristics of the data are clearly not the best match for the default hyperparameter combinations of XGBoost. This mismatch between the hyperparameters and the data characteristics of the specific problem affects the stability and prediction accuracy of the XGBoost ensemble algorithm. Therefore, to obtain an XGBoost ensemble model with better stability and higher prediction accuracy, it is necessary to optimize the hyperparameters. Therefore, the following section will optimize the XGBoost ensemble algorithm to improve the model's prediction accuracy.

[0177] Step 8: Optimize the hyperparameters of the XGBoost ensemble algorithm. Based on the data characteristics of parameter inversion, select the hyperparameters that can be tuned from the XGBoost ensemble algorithm. The hyperparameters that can be tuned are shown in Table 10.

[0178] Table 10

[0179]

[0180] First, the optimal range of hyperparameters is narrowed down using the method of controlling a single variable. The specific implementation process is as follows:

[0181] (1) For each hyperparameter, the hyperparameter values ​​are discretized to obtain a certain number of equal points;

[0182] (2) Substitute the discretized hyperparameter values ​​into the XGBoost ensemble algorithm model to obtain the coefficient of determination and generalization error of each model;

[0183] (3) Based on the variation law of the coefficient of determination and the generalization error, obtain the optimal range of values ​​for each hyperparameter.

[0184] Plot the coefficient of determination score curve and generalization error score curve for different hyperparameter values ​​in Table 10. Figure 9 The figure shows the coefficient of determination score curves for different values ​​of the hyperparameter n_estimators, as shown. Figure 10 The figure shows the generalization error score curves when the hyperparameter n_estimators takes different values.

[0185] Based on the score curves of the coefficient of determination and generalization error for the hyperparameter n_estimators at different values, when n_estimators is around 15, its coefficient of determination score reaches its peak, and its generalization error score also reaches its minimum value. As the hyperparameter continues to increase, both its coefficient of determination and generalization error remain unchanged. Therefore, based on the score curves of the coefficient of determination and generalization error, the optimization search range for this hyperparameter can be narrowed to [5, 20]. When the value of this hyperparameter is greater than 20, it will not only fail to increase the accuracy of the model and reduce the error, but will also cause overfitting, increase the complexity of the model, and reduce the running efficiency of the model.

[0186] The range of hyperparameter values ​​for the XGBoost ensemble algorithm, determined by the single-variable control method, is shown in Table 11.

[0187] Table 11

[0188]

[0189] Bayesian optimization is used to find the optimal combination of hyperparameters. Based on Bayes' theorem, Bayesian optimization consists of two main parts: a probabilistic surrogate function and a data acquisition function. To find the extrema of the objective function, the probabilistic surrogate function is updated after each data acquisition. Based on the feedback from the probabilistic surrogate function, the acquisition function acquires data in regions where extrema are likely to occur. This process is repeated iteratively after each acquisition. The Bayesian optimization algorithm makes full use of historical computational data, avoiding many detours, and its computational cost is much lower than other algorithms, making it a smarter optimization method.

[0190] The process of using Bayesian optimization to find the optimal hyperparameters is as follows:

[0191] (1) Establish hyperparameter samples as x = [x1, x2, x3, ..., x n The objective function of the XGBoost model takes the value y = [y1, y2, y3, ..., y]. n ], y follows a multidimensional normal distribution;

[0192] (2) Establish the kernel function for the Gaussian process:

[0193]

[0194] in, and l are the parameters of the kernel function;

[0195] (3) Based on the kernel function, establish the kernel vector:

[0196]

[0197] (4) Establish the likelihood function of y:

[0198]

[0199] Where μ is the mean of x;

[0200] The optimal kernel matrix is ​​obtained by minimizing the likelihood function using gradient descent and calculating the parameters of the kernel function.

[0201] (5) Based on the optimal kernel matrix, establish a Gaussian process between the posterior probability and prior probability of y and the input variables:

[0202]

[0203] Where y * K is the value of the output variable obtained in the next iteration. * This is the kernel vector for the next iteration;

[0204] The Gaussian process described above establishes the functional relationship between the objective function value and the hyperparameters of the model.

[0205] (6) Using the above Gaussian process as a probabilistic surrogate model, establish POI (Probability of Improvement) as the learning function:

[0206]

[0207] Where f(X) is the objective function value of X, i.e., the Gaussian process derived above; f(X) + ) represents the current optimal objective function value, i.e., the posterior distribution of f(X); μ(x) and σ(x) are the mean and variance of the objective function obtained by the Gaussian process; ξ is the trade-off coefficient, which controls whether the optimization direction of the hyperparameters is biased towards "development" or "search"; based on the above POI learned function, the Monte Carlo method is applied to search for the above optimal hyperparameters.

[0208] (7) Substitute the optimal hyperparameters searched based on the learned function into the probabilistic proxy model to establish a new Gaussian process;

[0209] (8) Using the new Gaussian process as the probabilistic surrogate model, establish a new POI and apply the Monte Carlo method to perform a new round of hyperparameter optimization.

[0210] (9) Repeat the above operation until the iteration stopping condition is met to obtain the final optimal hyperparameter combination.

[0211] Based on Python, the Bayesian Optimization algorithm was called from the bayes_opt module. The hyperparameter combinations output after 30 rounds of iterative optimization calculation are shown in Table 12.

[0212] Table 12

[0213]

[0214]

[0215] Step 9: Define the hyperparameters of the XGBoost ensemble algorithm using the hyperparameter combinations output in Table 12. Then, use 80% of the dataset in the database as the training set to train the XGBoost ensemble algorithm with defined hyperparameters, and use the remaining 20% ​​of the data as the test set to test the model. The Poisson ratio parameter inversion analysis evaluation index is shown in Table 13.

[0216] Table 13

[0217]

[0218] A bar chart comparing some predicted values ​​with actual values ​​is shown below. Figure 11 As shown.

[0219] To verify the practicality of the predicted parameters, some of the predicted parameters were substituted into the numerical simulation model. The displacements generated by the predicted parameters were compared with those generated by the actual parameters. An error of less than 10% proves the usability of the predicted parameters. The comparison of the displacements generated by the predicted parameters and the actual parameters is shown in Table 14. It can be seen that the errors of all examples are less than 10%, indicating that the tunnel surrounding rock parameter inversion model based on the optimized XGBoost ensemble algorithm has practical value.

[0220] Table 14

[0221]

[0222]

Claims

1. A method for intelligent inversion analysis of tunnel surrounding rock parameters based on the XGBoost optimization algorithm, characterized in that, Includes the following steps: S1: Establish a three-dimensional finite element numerical simulation model for tunnel excavation and support, and obtain a sample library for inverting surrounding rock parameters; S2: Perform feature correlation and sensitivity analysis on displacement and parameter data in the sample library, evaluate the feasibility of each parameter as the parameter to be inverted, and evaluate the rationality of displacement feature combination. S3: Use the CART algorithm to perform parameter inversion analysis and select displacement features; S4: Using the CART algorithm as the base learner, establish an XGBoost integrated algorithm model for intelligent inversion analysis of tunnel surrounding rock parameters; S5: Apply the single-variable control method and Bayesian optimization method to perform hyperparameter optimization, and obtain the XGBoost integrated algorithm model for intelligent inversion analysis of optimized tunnel surrounding rock parameters; S6: Input the displacement features into the trained XGBoost model for intelligent inversion analysis of surrounding rock parameters to obtain the predicted values ​​of surrounding rock parameters.

2. The intelligent inversion analysis method for tunnel surrounding rock parameters based on the XGBoost optimization algorithm according to claim 1, characterized in that, In S1, orthogonal experimental design is carried out on the surrounding rock parameters of the rock strata where the tunnel is located to obtain multiple sets of surrounding rock parameter combinations; multiple sets of surrounding rock parameters are input into the numerical simulation calculation model of tunnel excavation and support, and the displacement values ​​of the surrounding rock deformation stabilization stage are extracted to obtain a sample library for surrounding rock parameter inversion.

3. The intelligent inversion analysis method for tunnel surrounding rock parameters based on the XGBoost optimization algorithm according to claim 1, characterized in that, In S2, the feature correlation analysis uses the Pearson correlation coefficient for calculation. The Pearson correlation coefficient between two variables X and Y is the quotient of their covariance and standard deviation. The closer the absolute value of the Pearson correlation coefficient is to 1, the stronger the linear relationship between the two features. The formula for the Pearson correlation coefficient is: 。 4. The intelligent inversion analysis method for tunnel surrounding rock parameters based on the XGBoost optimization algorithm according to claim 1, characterized in that, In S2, the average value of the parameters of the rock strata where the tunnel is located, 90% of the average value of the parameters, and 110% of the average value of the parameters are taken and substituted into the numerical simulation calculation model to calculate the corresponding displacement. By comparing the displacement changes when different parameters change by the same proportion, the parameter sensitivity of the surrounding rock is judged.

5. The intelligent inversion analysis method for tunnel surrounding rock parameters based on the XGBoost optimization algorithm according to claim 1, characterized in that, In S3, the CART algorithm is a binary decision tree model, and the specific calculation process is as follows: (1) Select the segmentation variable j, which is the characteristic value of each displacement; (2) Select the dividing point s according to the dividing principle; (3) For any fractional point s of any fractional variable i The training sample x is split into two sample subsets R. 1i and R 2i ; in, , j=1~n; , j=1~n; i is the number of the split point; (4) Let c1 and c2 be R respectively 1i and R 2i The average value of the target value y in the middle sample: ; ; (5) Based on the principle of minimizing the squared error, the optimal splitting variables and splitting points are selected: (6) Based on the selected splitting variable j and splitting point s, the training samples are divided into two sample subsets, R1 and R2; (7) In the sample subsets R1 and R2, repeat steps (1) to (6) to further divide the subsets R1 and R2 into smaller subsets; (8) Repeat steps (1) to (7) to further divide the subset until the termination condition is met; (9) Finally, through the above division, the training sample input space is divided into R1, R2, R3, ... R m A production decision tree with m leaf nodes: 。 6. The intelligent inversion analysis method for tunnel surrounding rock parameters based on the XGBoost optimization algorithm according to claim 1, characterized in that, In S4, using the CART algorithm as the base learner, the XGBoost ensemble algorithm model is constructed as follows: (1) The CART algorithm is iterated. Each iteration fits the residual of the CART decision tree obtained in the previous iteration to minimize the residual. The formula for calculating the residual is: in, For residuals, The true value of the sample output is given by t, where t is the number of iterations. For the input variables of the sample, This represents the predicted value of the CART decision tree model at the (t-1)th iteration; (2) Perform the 0th iteration: in, Let L be the predicted value from the 0th iteration, and L be the objective function. Let c be the true value of the sample output variable, and c be the parameter of the objective function, which is the average value of the sample output variable. (3) Perform the first iteration: Use the difference between the true value of the sample output variable and the average value c of the sample output variable as the sample output value. Based on the specified objective function, the CART algorithm is applied to construct the first CART decision tree model, obtaining the predicted values ​​of the sample output variables. ; (4) Perform the second iteration: using the original output values ​​of the first iteration sample. Compared with the predicted values ​​of the output variables of the first iteration sample The difference is used as the output value of the sample. Based on the specified objective function, the CART algorithm is applied to construct a second CART decision tree model to obtain the predicted values ​​of the sample output variables. ; (5) Following the same method as the second iteration, construct the 3rd, 4th, ..., Mth CART decision tree models to obtain the predicted values ​​of the sample output variables. ,..., Continue until the residuals meet the requirements; (6) Based on the above M CART decision tree models, calculate the output value of the predicted sample according to the following formula: 。 7. The intelligent inversion analysis method for tunnel surrounding rock parameters based on the XGBoost optimization algorithm according to claim 6, characterized in that, In S4, the objective function of the XGBoost algorithm is obtained by adding a regularization term to the conventional squared error loss function: Where T is the number of leaf nodes; and These are the regularization coefficients for the leaf nodes and their weights, respectively. It is the sum of the first derivatives of all samples under leaf node j; This is the first derivative of the loss function with respect to the predicted values ​​of the sample output variables; It is the sum of the second derivatives of the j samples under the leaf node; This is the second derivative of the loss function with respect to the predicted values ​​of the sample output variables; The sample set under leaf node j; It is the prediction value given in the (t-1)th round after the accumulation of the learner in the first t-1 rounds; is the original value of the sample output variable; i is the sample number.

8. The intelligent inversion analysis method for tunnel surrounding rock parameters based on the XGBoost optimization algorithm according to claim 1, characterized in that, In S5, the range of hyperparameter values ​​is narrowed down using the single-variable control method. Within this narrowed range, the optimal values ​​of the hyperparameters are further determined using Bayesian optimization. The process of optimizing the hyperparameters using Bayesian optimization is as follows: (1) Establish hyperparameter samples as x=[x1, x2, x3,...,x n The objective function of the XGBoost ensemble algorithm model takes the value y = [y1, y2, y3, ..., y]. n ], y follows a multidimensional normal distribution; (2) Establish the kernel function of the Gaussian process: in, and These are the parameters of the kernel function; (3) Based on the kernel function, establish the kernel vector: (4) Establish the likelihood function of y: in, Let x be the mean of x; The optimal kernel matrix is ​​obtained by minimizing the likelihood function using gradient descent and calculating the parameters of the kernel function. (5) Based on the optimal kernel matrix, establish a Gaussian process between the posterior probability and prior probability of y and the input variables: [ in To assign values ​​to the output variables obtained in the next iteration, This is the kernel vector for the next iteration; The above Gaussian process establishes the functional relationship between the objective function value and hyperparameters of the XGBoost ensemble algorithm model; (6) Using the above Gaussian process as a probabilistic surrogate model, establish POI as the learned function: in, for The objective function value, i.e., the Gaussian process derived above; The current optimal objective function value is, The posterior distribution of; and Let be the mean and variance of the objective function obtained from the Gaussian process; The trade-off coefficient controls whether the optimization direction of the hyperparameters leans towards development or search; based on the POI learned function, the Monte Carlo method is applied to search for the optimal hyperparameters; (7) Substitute the optimal hyperparameters obtained from the learned function search into the probabilistic surrogate model to establish a new Gaussian process; (8) Using the new Gaussian process as the probabilistic surrogate model, establish a new POI and apply the Monte Carlo method to perform a new round of hyperparameter optimization; Repeat the above steps until the iteration stopping condition is met to obtain the final optimal combination of hyperparameters.

Images